Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+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
1463 void
1465};
1466
1467
1468
1486template <int n_components_,
1487 int dim,
1488 int spacedim = dim,
1489 typename Number = double>
1491 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1492{
1493public:
1494 static constexpr unsigned int dimension = dim;
1495 static constexpr unsigned int n_components = n_components_;
1496
1497 using number_type = Number;
1498
1501 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1503 using ETT = typename internal::FEPointEvaluation::
1504 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1505 using value_type = typename ETT::value_type;
1506 using scalar_value_type = typename ETT::scalar_value_type;
1507 using vectorized_value_type = typename ETT::vectorized_value_type;
1508 using unit_gradient_type = typename ETT::unit_gradient_type;
1509 using gradient_type = typename ETT::real_gradient_type;
1511 typename ETT::interface_vectorized_unit_gradient_type;
1512
1519 const bool is_interior = true,
1520 const unsigned int first_selected_component = 0);
1521
1526 void
1527 reinit(const unsigned int cell_index, const unsigned int face_number);
1528
1533 void
1534 reinit(const unsigned int face_index);
1535
1547 template <std::size_t stride_view>
1548 void
1549 evaluate(
1551 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1552
1564 void
1565 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1566 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1567
1590 template <std::size_t stride_view>
1591 void
1593 const EvaluationFlags::EvaluationFlags &integration_flags,
1594 const bool sum_into_values = false);
1595
1618 void
1619 integrate(const ArrayView<ScalarNumber> &solution_values,
1620 const EvaluationFlags::EvaluationFlags &integration_flags,
1621 const bool sum_into_values = false);
1622
1645 template <std::size_t stride_view>
1646 void
1648 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1649 const EvaluationFlags::EvaluationFlags &integration_flags,
1650 const bool sum_into_values = false);
1651
1674 void
1675 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1676 const EvaluationFlags::EvaluationFlags &integration_flags,
1677 const bool sum_into_values = false);
1678
1685 template <int stride_face_dof = VectorizedArrayType::size()>
1686 void
1687 evaluate_in_face(const ScalarNumber *face_dof_values,
1688 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1689
1696 template <int stride_face_dof = VectorizedArrayType::size()>
1697 void
1698 integrate_in_face(ScalarNumber *face_dof_values,
1699 const EvaluationFlags::EvaluationFlags &integration_flags,
1700 const bool sum_into_values = false);
1701
1708 normal_vector(const unsigned int point_index) const;
1709
1710private:
1711 static constexpr std::size_t n_lanes_user_interface =
1713 static constexpr std::size_t n_lanes_internal =
1715 static constexpr std::size_t stride =
1717
1718 template <bool is_linear, std::size_t stride_view>
1719 void
1722 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1723
1724 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1725 void
1727 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1728 const EvaluationFlags::EvaluationFlags &integration_flags,
1729 const bool sum_into_values);
1730
1735 template <bool is_linear, int stride_face_dof>
1736 void
1737 do_evaluate_in_face(const ScalarNumber *face_dof_values,
1738 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1739
1744 template <bool do_JxW, bool is_linear, int stride_face_dof>
1745 void
1747 ScalarNumber *face_dof_values,
1748 const EvaluationFlags::EvaluationFlags &integration_flags,
1749 const bool sum_into_values);
1750};
1751
1752
1753
1754// ----------------------- template and inline function ----------------------
1755
1756
1757template <int n_components_, int dim, int spacedim, typename Number>
1761 const UpdateFlags update_flags,
1762 const unsigned int first_selected_component)
1763 : n_q_batches(numbers::invalid_unsigned_int)
1764 , n_q_points(numbers::invalid_unsigned_int)
1765 , n_q_points_scalar(numbers::invalid_unsigned_int)
1766 , mapping(&mapping)
1767 , fe(&fe)
1768 , JxW_ptr(nullptr)
1769 , update_flags(update_flags)
1770 , mapping_info_on_the_fly(
1771 std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1772 mapping,
1773 update_flags))
1774 , mapping_info(mapping_info_on_the_fly.get())
1775 , current_cell_index(numbers::invalid_unsigned_int)
1776 , current_face_number(numbers::invalid_unsigned_int)
1777 , must_reinitialize_pointers(false)
1778 , is_interior(true)
1779{
1780 setup(first_selected_component);
1781}
1782
1783
1784
1785template <int n_components_, int dim, int spacedim, typename Number>
1790 const unsigned int first_selected_component,
1791 const bool is_interior)
1792 : n_q_batches(numbers::invalid_unsigned_int)
1793 , n_q_points(numbers::invalid_unsigned_int)
1794 , n_q_points_scalar(numbers::invalid_unsigned_int)
1795 , mapping(&mapping_info.get_mapping())
1796 , fe(&fe)
1797 , JxW_ptr(nullptr)
1798 , update_flags(mapping_info.get_update_flags())
1799 , mapping_info(&mapping_info)
1800 , current_cell_index(numbers::invalid_unsigned_int)
1801 , current_face_number(numbers::invalid_unsigned_int)
1802 , must_reinitialize_pointers(true)
1803 , is_interior(is_interior)
1804{
1805 setup(first_selected_component);
1806}
1807
1808
1809
1810template <int n_components_, int dim, int spacedim, typename Number>
1814 : n_q_batches(other.n_q_batches)
1815 , n_q_points(other.n_q_points)
1816 , n_q_points_scalar(other.n_q_points_scalar)
1817 , mapping(other.mapping)
1818 , fe(other.fe)
1819 , poly(other.poly)
1820 , use_linear_path(other.use_linear_path)
1821 , renumber(other.renumber)
1822 , solution_renumbered(other.solution_renumbered)
1823 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1824 , values(other.values)
1825 , gradients(other.gradients)
1826 , dofs_per_component(other.dofs_per_component)
1827 , dofs_per_component_face(other.dofs_per_component_face)
1828 , component_in_base_element(other.component_in_base_element)
1829 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1830 , update_flags(other.update_flags)
1831 , fe_values(other.fe_values)
1832 , mapping_info_on_the_fly(
1833 other.mapping_info_on_the_fly ?
1835 *mapping,
1836 update_flags) :
1837 nullptr)
1838 , mapping_info(other.mapping_info)
1839 , current_cell_index(other.current_cell_index)
1840 , current_face_number(other.current_face_number)
1841 , fast_path(other.fast_path)
1842 , must_reinitialize_pointers(true)
1843 , is_interior(other.is_interior)
1844{}
1845
1846
1847
1848template <int n_components_, int dim, int spacedim, typename Number>
1851 : n_q_batches(other.n_q_batches)
1852 , n_q_points(other.n_q_points)
1853 , n_q_points_scalar(other.n_q_points_scalar)
1854 , mapping(other.mapping)
1855 , fe(other.fe)
1856 , poly(other.poly)
1857 , use_linear_path(other.use_linear_path)
1858 , renumber(other.renumber)
1859 , solution_renumbered(other.solution_renumbered)
1860 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1861 , values(other.values)
1862 , gradients(other.gradients)
1863 , dofs_per_component(other.dofs_per_component)
1864 , dofs_per_component_face(other.dofs_per_component_face)
1865 , component_in_base_element(other.component_in_base_element)
1866 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1867 , update_flags(other.update_flags)
1868 , fe_values(other.fe_values)
1869 , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1870 , mapping_info(other.mapping_info)
1871 , current_cell_index(other.current_cell_index)
1872 , current_face_number(other.current_face_number)
1873 , fast_path(other.fast_path)
1874 , must_reinitialize_pointers(other.must_reinitialize_pointers)
1875 , is_interior(other.is_interior)
1876{}
1877
1878
1879
1880template <int n_components_, int dim, int spacedim, typename Number>
1881void
1883 const unsigned int first_selected_component)
1884{
1885 AssertIndexRange(first_selected_component + n_components,
1886 fe->n_components() + 1);
1887
1888 shapes.reserve(100);
1889
1890 bool same_base_element = true;
1891 unsigned int base_element_number = 0;
1892 component_in_base_element = 0;
1893 unsigned int component = 0;
1894 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1895 if (component + fe->element_multiplicity(base_element_number) >
1896 first_selected_component)
1897 {
1898 if (first_selected_component + n_components >
1899 component + fe->element_multiplicity(base_element_number))
1900 same_base_element = false;
1901 component_in_base_element = first_selected_component - component;
1902 break;
1903 }
1904 else
1905 component += fe->element_multiplicity(base_element_number);
1906
1909 *fe, base_element_number) &&
1910 same_base_element)
1911 {
1912 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1913 renumber = shape_info.lexicographic_numbering;
1914 dofs_per_component = shape_info.dofs_per_component_on_cell;
1915 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1917 fe->base_element(base_element_number));
1918
1919 bool is_lexicographic = true;
1920 for (unsigned int i = 0; i < renumber.size(); ++i)
1921 if (i != renumber[i])
1922 is_lexicographic = false;
1923
1924 if (is_lexicographic)
1925 renumber.clear();
1926
1927 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1928 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1929 poly[1].value(1.) == 1.) &&
1930 (fe->n_components() == n_components);
1931
1932 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1933 const unsigned int size_cell = dofs_per_component * n_components;
1934 scratch_data_scalar.resize(size_face + size_cell);
1935
1936 solution_renumbered.resize(dofs_per_component);
1937 solution_renumbered_vectorized.resize(dofs_per_component);
1938
1939 fast_path = true;
1940 }
1941 else
1942 {
1943 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1944 for (unsigned int d = 0; d < n_components; ++d)
1945 {
1946 const unsigned int component = first_selected_component + d;
1947 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1948 {
1949 const bool is_primitive =
1950 fe->is_primitive() || fe->is_primitive(i);
1951 if (is_primitive)
1952 nonzero_shape_function_component[i][d] =
1953 (component == fe->system_to_component_index(i).first);
1954 else
1955 nonzero_shape_function_component[i][d] =
1956 (fe->get_nonzero_components(i)[component] == true);
1957 }
1958 }
1959
1960 fast_path = false;
1961 }
1962}
1963
1964
1965
1966template <int n_components_, int dim, int spacedim, typename Number>
1967template <bool is_face, bool is_linear>
1968inline void
1970{
1971 const unsigned int geometry_index =
1972 mapping_info->template compute_geometry_index_offset<is_face>(
1973 current_cell_index, current_face_number);
1974
1975 cell_type = mapping_info->get_cell_type(geometry_index);
1976
1977 const_cast<unsigned int &>(n_q_points_scalar) =
1978 mapping_info->get_n_q_points_unvectorized(geometry_index);
1979
1980 // round up n_q_points_scalar / n_lanes_internal
1981 const_cast<unsigned int &>(n_q_batches) =
1982 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
1983
1984 const unsigned int n_q_points_before = n_q_points;
1985
1986 const_cast<unsigned int &>(n_q_points) =
1987 (stride == 1) ? n_q_batches : n_q_points_scalar;
1988
1989 if (n_q_points != n_q_points_before)
1990 {
1991 if (update_flags & update_values)
1992 values.resize(n_q_points);
1993 if (update_flags & update_gradients)
1994 gradients.resize(n_q_points);
1995 }
1996
1997 if (n_q_points == 0)
1998 return;
1999
2000 // set unit point pointer
2001 const unsigned int unit_point_offset =
2002 mapping_info->compute_unit_point_index_offset(geometry_index);
2003
2004 if (is_face)
2005 unit_point_faces_ptr =
2006 mapping_info->get_unit_point_faces(unit_point_offset);
2007 else
2008 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2009
2010 // set data pointers
2011 const unsigned int data_offset =
2012 mapping_info->compute_data_index_offset(geometry_index);
2013 const unsigned int compressed_data_offset =
2014 mapping_info->compute_compressed_data_index_offset(geometry_index);
2015#ifdef DEBUG
2016 const UpdateFlags update_flags_mapping =
2017 mapping_info->get_update_flags_mapping();
2018 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2019 real_point_ptr = mapping_info->get_real_point(data_offset);
2020 if (update_flags_mapping & UpdateFlags::update_jacobians)
2021 jacobian_ptr =
2022 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2023 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2024 inverse_jacobian_ptr =
2025 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2026 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2027 normal_ptr = mapping_info->get_normal_vector(data_offset);
2028 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2029 JxW_ptr = mapping_info->get_JxW(data_offset);
2030#else
2031 real_point_ptr = mapping_info->get_real_point(data_offset);
2032 jacobian_ptr =
2033 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2034 inverse_jacobian_ptr =
2035 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2036 normal_ptr = mapping_info->get_normal_vector(data_offset);
2037 JxW_ptr = mapping_info->get_JxW(data_offset);
2038#endif
2039
2040 if (!is_linear && fast_path)
2041 {
2042 const std::size_t n_shapes = poly.size();
2043 if (is_face)
2044 shapes_faces.resize_fast(n_q_batches * n_shapes);
2045 else
2046 shapes.resize_fast(n_q_batches * n_shapes);
2047
2048 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2049 if (is_face)
2050 {
2051 if (dim > 1)
2052 {
2054 shapes_faces.data() + qb * n_shapes,
2055 poly,
2056 unit_point_faces_ptr[qb],
2057 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2058 }
2059 }
2060 else
2061 {
2062 if (update_flags & UpdateFlags::update_gradients)
2063 {
2064 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2065 poly,
2066 unit_point_ptr[qb],
2067 1);
2068 }
2069 else if (qb + 1 < n_q_batches)
2070 {
2071 // Use function with reduced overhead to compute for two
2072 // points at once
2074 shapes.data() + qb * n_shapes,
2075 poly,
2076 unit_point_ptr[qb],
2077 unit_point_ptr[qb + 1]);
2078 ++qb;
2079 }
2080 else
2081 {
2082 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2083 poly,
2084 unit_point_ptr[qb],
2085 0);
2086 }
2087 }
2088 }
2089}
2090
2091
2092
2093template <int n_components_, int dim, int spacedim, typename Number>
2094inline const typename FEPointEvaluationBase<n_components_,
2095 dim,
2096 spacedim,
2097 Number>::value_type &
2099 const unsigned int point_index) const
2100{
2101 AssertIndexRange(point_index, values.size());
2102 return values[point_index];
2103}
2104
2105
2106
2107template <int n_components_, int dim, int spacedim, typename Number>
2108inline const typename FEPointEvaluationBase<n_components_,
2109 dim,
2110 spacedim,
2111 Number>::gradient_type &
2113 const unsigned int point_index) const
2114{
2115 AssertIndexRange(point_index, gradients.size());
2116 return gradients[point_index];
2117}
2118
2119
2120
2121template <int n_components_, int dim, int spacedim, typename Number>
2122inline void
2124 const value_type &value,
2125 const unsigned int point_index)
2126{
2127 AssertIndexRange(point_index, n_q_points);
2128 values[point_index] = value;
2129}
2130
2131
2132
2133template <int n_components_, int dim, int spacedim, typename Number>
2134inline void
2136 const gradient_type &gradient,
2137 const unsigned int point_index)
2138{
2139 AssertIndexRange(point_index, n_q_points);
2140 gradients[point_index] = gradient;
2141}
2142
2143
2144
2145template <int n_components_, int dim, int spacedim, typename Number>
2148 const unsigned int point_index) const
2149{
2150 AssertIndexRange(point_index, n_q_points);
2151 Assert(jacobian_ptr != nullptr,
2153 ExcFEPointEvaluationAccessToUninitializedMappingField(
2154 "update_jacobians"));
2155 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2156 GeometryType::affine ?
2157 0 :
2158 point_index];
2159}
2160
2161
2162
2163template <int n_components_, int dim, int spacedim, typename Number>
2166 const unsigned int point_index) const
2167{
2168 AssertIndexRange(point_index, n_q_points);
2169 Assert(inverse_jacobian_ptr != nullptr,
2171 ExcFEPointEvaluationAccessToUninitializedMappingField(
2172 "update_inverse_jacobians"));
2173 return inverse_jacobian_ptr
2174 [cell_type <=
2176 0 :
2177 point_index];
2178}
2179
2180
2181
2182template <int n_components_, int dim, int spacedim, typename Number>
2183inline Number
2185 const unsigned int point_index) const
2186{
2187 AssertIndexRange(point_index, n_q_points);
2188 Assert(JxW_ptr != nullptr,
2190 ExcFEPointEvaluationAccessToUninitializedMappingField(
2191 "update_JxW_values"));
2192 return JxW_ptr[point_index];
2193}
2194
2195
2196
2197template <int n_components_, int dim, int spacedim, typename Number>
2200 const unsigned int point_index) const
2201{
2202 return quadrature_point(point_index);
2203 ;
2204}
2205
2206
2207
2208template <int n_components_, int dim, int spacedim, typename Number>
2211 const unsigned int point_index) const
2212{
2213 AssertIndexRange(point_index, n_q_points);
2214 Assert(real_point_ptr != nullptr,
2216 ExcFEPointEvaluationAccessToUninitializedMappingField(
2217 "update_quadrature_points"));
2218 return real_point_ptr[point_index];
2219}
2220
2221
2222
2223template <int n_components_, int dim, int spacedim, typename Number>
2224inline Point<dim, Number>
2226 const unsigned int point_index) const
2227{
2228 AssertIndexRange(point_index, n_q_points);
2229 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2230 Point<dim, Number> unit_point;
2231 for (unsigned int d = 0; d < dim; ++d)
2233 unit_point_ptr[point_index / stride][d], point_index % stride);
2234 return unit_point;
2235}
2236
2237
2238
2239template <int n_components_, int dim, int spacedim, typename Number>
2246
2247
2248
2249template <int n_components_, int dim, int spacedim, typename Number>
2253 const unsigned int first_selected_component)
2254 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2255 mapping_info,
2256 fe,
2257 first_selected_component)
2258{}
2259
2260
2261
2262template <int n_components_, int dim, int spacedim, typename Number>
2264 const Mapping<dim, spacedim> &mapping,
2266 const UpdateFlags update_flags,
2267 const unsigned int first_selected_component)
2268 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2269 mapping,
2270 fe,
2271 update_flags,
2272 first_selected_component)
2273{}
2274
2275
2276
2277template <int n_components_, int dim, int spacedim, typename Number>
2278inline void
2281{
2282 this->current_cell_index = numbers::invalid_unsigned_int;
2283 this->current_face_number = numbers::invalid_unsigned_int;
2284
2285 if (this->use_linear_path)
2286 this->template do_reinit<false, true>();
2287 else
2288 this->template do_reinit<false, false>();
2289}
2290
2291
2292
2293template <int n_components_, int dim, int spacedim, typename Number>
2294inline void
2296{
2297 internal_reinit_single_cell_state_mapping_info();
2298 this->must_reinitialize_pointers = false;
2299}
2300
2301
2302
2303template <int n_components_, int dim, int spacedim, typename Number>
2304inline void
2307 const ArrayView<const Point<dim>> &unit_points)
2308{
2309 // reinit is only allowed for mapping computation on the fly
2310 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2312
2313 this->mapping_info->reinit(cell, unit_points);
2314 this->must_reinitialize_pointers = false;
2315
2316 if (!this->fast_path)
2317 {
2318 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2319 *this->mapping,
2320 *this->fe,
2322 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2323 this->update_flags);
2324 this->fe_values->reinit(cell);
2325 }
2326
2327 if (this->use_linear_path)
2328 this->template do_reinit<false, true>();
2329 else
2330 this->template do_reinit<false, false>();
2331}
2332
2333
2334
2335template <int n_components_, int dim, int spacedim, typename Number>
2336inline void
2338 const unsigned int cell_index)
2339{
2340 this->current_cell_index = cell_index;
2341 this->current_face_number = numbers::invalid_unsigned_int;
2342 this->must_reinitialize_pointers = false;
2343
2344 if (this->use_linear_path)
2345 this->template do_reinit<false, true>();
2346 else
2347 this->template do_reinit<false, false>();
2348
2349 if (!this->fast_path)
2350 {
2351 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2352
2353 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2354 for (unsigned int d = 0; d < dim; ++d)
2355 unit_points[v][d] =
2356 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2357
2358 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2359 *this->mapping,
2360 *this->fe,
2362 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2363 this->update_flags);
2364
2365 this->fe_values->reinit(
2366 this->mapping_info->get_cell_iterator(this->current_cell_index));
2367 }
2368}
2369
2370
2371
2372template <int n_components_, int dim, int spacedim, typename Number>
2373template <std::size_t stride_view>
2374void
2377 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2378{
2379 if (this->must_reinitialize_pointers)
2380 internal_reinit_single_cell_state_mapping_info();
2381
2382 if (this->n_q_points == 0)
2383 return;
2384
2385 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2386
2387 if (!((evaluation_flags & EvaluationFlags::values) ||
2388 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2389 return;
2390
2391 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2392 if (this->fast_path)
2393 {
2394 if (this->use_linear_path)
2395 evaluate_fast<true>(solution_values, evaluation_flags);
2396 else
2397 evaluate_fast<false>(solution_values, evaluation_flags);
2398 }
2399 else
2400 evaluate_slow(solution_values, evaluation_flags);
2401}
2402
2403
2404
2405template <int n_components_, int dim, int spacedim, typename Number>
2406void
2408 const ArrayView<const ScalarNumber> &solution_values,
2409 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2410{
2411 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2412 solution_values.size()),
2413 evaluation_flags);
2414}
2415
2416
2417
2418template <int n_components_, int dim, int spacedim, typename Number>
2419template <std::size_t stride_view>
2420void
2422 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2423 const EvaluationFlags::EvaluationFlags &integration_flags,
2424 const bool sum_into_values)
2425{
2426 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2427}
2428
2429
2430
2431template <int n_components_, int dim, int spacedim, typename Number>
2432void
2434 const ArrayView<ScalarNumber> &solution_values,
2435 const EvaluationFlags::EvaluationFlags &integration_flags,
2436 const bool sum_into_values)
2437{
2438 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2439 solution_values.size()),
2440 integration_flags,
2441 sum_into_values);
2442}
2443
2444
2445
2446template <int n_components_, int dim, int spacedim, typename Number>
2448 scalar_value_type
2450 const
2451{
2452 value_type return_value = {};
2453
2454 for (const auto point_index : this->quadrature_point_indices())
2455 return_value += values[point_index] * this->JxW(point_index);
2456
2457 return ETT::sum_value(return_value);
2458}
2459
2460
2461
2462template <int n_components_, int dim, int spacedim, typename Number>
2463unsigned int
2466{
2467 Assert(stride == 1,
2468 ExcMessage(
2469 "Calling this function only makes sense in fully vectorized mode."));
2470 if (q == n_q_batches - 1)
2471 {
2472 const unsigned int n_filled_lanes =
2473 n_q_points_scalar & (n_lanes_user_interface - 1);
2474 if (n_filled_lanes == 0)
2475 return n_lanes_user_interface;
2476 else
2477 return n_filled_lanes;
2478 }
2479 else
2480 return n_lanes_user_interface;
2481}
2482
2483
2484
2485template <int n_components_, int dim, int spacedim, typename Number>
2486template <std::size_t stride_view>
2487void
2489 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2490 const EvaluationFlags::EvaluationFlags &integration_flags,
2491 const bool sum_into_values)
2492{
2493 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2494}
2495
2496
2497
2498template <int n_components_, int dim, int spacedim, typename Number>
2499void
2501 const ArrayView<ScalarNumber> &solution_values,
2502 const EvaluationFlags::EvaluationFlags &integration_flags,
2503 const bool sum_into_values)
2504{
2505 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2506 solution_values.size()),
2507 integration_flags,
2508 sum_into_values);
2509}
2510
2511
2512
2513template <int n_components_, int dim, int spacedim, typename Number>
2514template <bool is_linear, std::size_t stride_view>
2515inline void
2518{
2519 const unsigned int dofs_per_comp =
2520 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2521
2522 for (unsigned int comp = 0; comp < n_components; ++comp)
2523 {
2524 const std::size_t offset =
2525 (this->component_in_base_element + comp) * dofs_per_comp;
2526
2527 if ((is_linear && n_components == 1) || this->renumber.empty())
2528 {
2529 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2530 ETT::read_value(solution_values[i + offset],
2531 comp,
2532 this->solution_renumbered[i]);
2533 }
2534 else
2535 {
2536 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2537 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2538 ETT::read_value(solution_values[renumber_ptr[i]],
2539 comp,
2540 this->solution_renumbered[i]);
2541 }
2542 }
2543}
2544
2545
2546
2547template <int n_components_, int dim, int spacedim, typename Number>
2548template <bool is_linear, std::size_t stride_view>
2549inline void
2552 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2553 const unsigned int n_shapes,
2554 const unsigned int qb,
2555 vectorized_value_type &value,
2557{
2558 if (evaluation_flags & EvaluationFlags::gradients)
2559 {
2560 std::array<vectorized_value_type, dim + 1> result;
2561 if constexpr (is_linear)
2562 {
2563 if constexpr (n_components == 1)
2564 result =
2566 dim,
2569 1,
2570 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2571 else
2572 result =
2574 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2575 }
2576 else
2578 dim,
2581 1,
2582 false>(this->shapes.data() + qb * n_shapes,
2583 n_shapes,
2584 this->solution_renumbered.data());
2585 gradient[0] = result[0];
2586 if (dim > 1)
2587 gradient[1] = result[1];
2588 if (dim > 2)
2589 gradient[2] = result[2];
2590 value = result[dim];
2591 }
2592 else
2593 {
2594 if constexpr (is_linear)
2595 {
2596 if constexpr (n_components == 1)
2598 dim,
2601 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2602 else
2604 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2605 }
2606 else
2607 value =
2611 false>(
2612 this->shapes.data() + qb * n_shapes,
2613 n_shapes,
2614 this->solution_renumbered.data());
2615 }
2616}
2617
2618
2619
2620template <int n_components_, int dim, int spacedim, typename Number>
2621template <bool is_linear, std::size_t stride_view>
2622inline void
2625 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2626{
2627 if (!(is_linear && n_components == 1))
2628 prepare_evaluate_fast<is_linear>(solution_values);
2629
2630 // loop over quadrature batches qb
2631 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2632
2633 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2634 {
2637
2638 compute_evaluate_fast<is_linear>(
2639 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2640
2641 if (evaluation_flags & EvaluationFlags::values)
2642 {
2643 for (unsigned int v = 0, offset = qb * stride;
2644 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2645 ++v, ++offset)
2646 ETT::set_value(value, v, this->values[offset]);
2647 }
2648 if (evaluation_flags & EvaluationFlags::gradients)
2649 {
2650 Assert(this->update_flags & update_gradients ||
2651 this->update_flags & update_inverse_jacobians,
2653
2654 for (unsigned int v = 0, offset = qb * stride;
2655 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2656 ++v, ++offset)
2657 {
2658 unit_gradient_type unit_gradient;
2659 ETT::set_gradient(gradient, v, unit_gradient);
2660 this->gradients[offset] =
2661 this->cell_type <=
2664 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2666 this
2667 ->inverse_jacobian_ptr[this->cell_type <=
2669 GeometryType::affine ?
2670 0 :
2671 offset]
2672 .transpose(),
2673 unit_gradient);
2674 }
2675 }
2676 }
2677}
2678
2679
2680
2681template <int n_components_, int dim, int spacedim, typename Number>
2682template <std::size_t stride_view>
2683inline void
2686 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2687{
2688 // slow path with FEValues
2689 Assert(this->fe_values.get() != nullptr,
2690 ExcMessage(
2691 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2692
2693 const std::size_t n_points = this->fe_values->get_quadrature().size();
2694
2695 if (evaluation_flags & EvaluationFlags::values)
2696 {
2697 this->values.resize(this->n_q_points);
2698 std::fill(this->values.begin(), this->values.end(), value_type());
2699 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2700 {
2701 const ScalarNumber value = solution_values[i];
2702 for (unsigned int d = 0; d < n_components; ++d)
2703 if (this->nonzero_shape_function_component[i][d] &&
2704 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2705 for (unsigned int qb = 0, q = 0; q < n_points;
2706 ++qb, q += n_lanes_user_interface)
2707 for (unsigned int v = 0;
2708 v < n_lanes_user_interface && q + v < n_points;
2709 ++v)
2710 ETT::access(this->values[qb],
2711 v,
2712 d,
2713 this->fe_values->shape_value(i, q + v) * value);
2714 else if (this->nonzero_shape_function_component[i][d])
2715 for (unsigned int qb = 0, q = 0; q < n_points;
2716 ++qb, q += n_lanes_user_interface)
2717 for (unsigned int v = 0;
2718 v < n_lanes_user_interface && q + v < n_points;
2719 ++v)
2720 ETT::access(this->values[qb],
2721 v,
2722 d,
2723 this->fe_values->shape_value_component(i,
2724 q + v,
2725 d) *
2726 value);
2727 }
2728 }
2729
2730 if (evaluation_flags & EvaluationFlags::gradients)
2731 {
2732 this->gradients.resize(this->n_q_points);
2733 std::fill(this->gradients.begin(),
2734 this->gradients.end(),
2735 gradient_type());
2736 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2737 {
2738 const ScalarNumber value = solution_values[i];
2739 for (unsigned int d = 0; d < n_components; ++d)
2740 if (this->nonzero_shape_function_component[i][d] &&
2741 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2742 for (unsigned int qb = 0, q = 0; q < n_points;
2743 ++qb, q += n_lanes_user_interface)
2744 for (unsigned int v = 0;
2745 v < n_lanes_user_interface && q + v < n_points;
2746 ++v)
2747 ETT::access(this->gradients[qb],
2748 v,
2749 d,
2750 this->fe_values->shape_grad(i, q + v) * value);
2751 else if (this->nonzero_shape_function_component[i][d])
2752 for (unsigned int qb = 0, q = 0; q < n_points;
2753 ++qb, q += n_lanes_user_interface)
2754 for (unsigned int v = 0;
2755 v < n_lanes_user_interface && q + v < n_points;
2756 ++v)
2757 ETT::access(
2758 this->gradients[qb],
2759 v,
2760 d,
2761 this->fe_values->shape_grad_component(i, q + v, d) * value);
2762 }
2763 }
2764}
2765
2766
2767
2768template <int n_components_, int dim, int spacedim, typename Number>
2769template <bool is_linear>
2770inline void
2772 const EvaluationFlags::EvaluationFlags &integration_flags,
2773 const unsigned int n_shapes,
2774 const unsigned int qb,
2775 const vectorized_value_type value,
2777 vectorized_value_type *solution_values_vectorized_linear)
2778{
2779 if (integration_flags & EvaluationFlags::gradients)
2781 is_linear,
2782 dim,
2784 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2785 n_shapes,
2786 &value,
2787 gradient,
2788 is_linear ?
2789 solution_values_vectorized_linear :
2790 this->solution_renumbered_vectorized.data(),
2791 this->unit_point_ptr[qb],
2792 qb != 0);
2793 else
2795 dim,
2798 this->shapes.data() + qb * n_shapes,
2799 n_shapes,
2800 value,
2801 is_linear ? solution_values_vectorized_linear :
2802 this->solution_renumbered_vectorized.data(),
2803 this->unit_point_ptr[qb],
2804 qb != 0);
2805}
2806
2807
2808
2809template <int n_components_, int dim, int spacedim, typename Number>
2810template <bool is_linear, std::size_t stride_view>
2811inline void
2813 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2814 vectorized_value_type *solution_values_vectorized_linear,
2815 const bool sum_into_values)
2816{
2817 if (!sum_into_values && this->fe->n_components() > n_components)
2818 for (unsigned int i = 0; i < solution_values.size(); ++i)
2819 solution_values[i] = 0;
2820
2821 const unsigned int dofs_per_comp =
2822 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2823
2824 for (unsigned int comp = 0; comp < n_components; ++comp)
2825 {
2826 const std::size_t offset =
2827 (this->component_in_base_element + comp) * dofs_per_comp;
2828
2829 if (is_linear || this->renumber.empty())
2830 {
2831 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2832 if (sum_into_values)
2833 solution_values[i + offset] +=
2834 ETT::sum_value(comp,
2835 is_linear ?
2836 *(solution_values_vectorized_linear + i) :
2837 this->solution_renumbered_vectorized[i]);
2838 else
2839 solution_values[i + offset] =
2840 ETT::sum_value(comp,
2841 is_linear ?
2842 *(solution_values_vectorized_linear + i) :
2843 this->solution_renumbered_vectorized[i]);
2844 }
2845 else
2846 {
2847 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2848 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2849 if (sum_into_values)
2850 solution_values[renumber_ptr[i]] +=
2851 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2852 else
2853 solution_values[renumber_ptr[i]] =
2854 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2855 }
2856 }
2857}
2858
2859
2860
2861template <int n_components_, int dim, int spacedim, typename Number>
2862template <bool do_JxW, bool is_linear, std::size_t stride_view>
2863inline void
2865 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2866 const EvaluationFlags::EvaluationFlags &integration_flags,
2867 const bool sum_into_values)
2868{
2869 // zero out lanes of incomplete last quadrature point batch
2870 if constexpr (stride == 1)
2871 if (const unsigned int n_filled_lanes =
2872 this->n_q_points_scalar & (n_lanes_internal - 1);
2873 n_filled_lanes > 0)
2874 {
2875 if (integration_flags & EvaluationFlags::values)
2876 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2877 ETT::set_zero_value(this->values.back(), v);
2878 if (integration_flags & EvaluationFlags::gradients)
2879 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2880 ETT::set_zero_gradient(this->gradients.back(), v);
2881 }
2882
2883 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
2884 solution_values_vectorized_linear = {};
2885
2886 // loop over quadrature batches qb
2887 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2888
2889 const bool cartesian_cell =
2891 const bool affine_cell =
2893 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2894 {
2895 vectorized_value_type value = {};
2897
2898 if (integration_flags & EvaluationFlags::values)
2899 for (unsigned int v = 0, offset = qb * stride;
2900 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2901 ++v, ++offset)
2902 ETT::get_value(value,
2903 v,
2904 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
2905 this->values[offset]);
2906
2907 if (integration_flags & EvaluationFlags::gradients)
2908 for (unsigned int v = 0, offset = qb * stride;
2909 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2910 ++v, ++offset)
2911 {
2912 const gradient_type grad_w =
2913 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
2914 this->gradients[offset];
2915 ETT::get_gradient(
2916 gradient,
2917 v,
2918 cartesian_cell ?
2919 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
2920 grad_w) :
2922 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
2923 grad_w));
2924 }
2925
2926 compute_integrate_fast<is_linear>(
2927 integration_flags,
2928 n_shapes,
2929 qb,
2930 value,
2931 gradient,
2932 solution_values_vectorized_linear.data());
2933 }
2934
2935 // add between the lanes and write into the result
2936 finish_integrate_fast<is_linear>(solution_values,
2937 solution_values_vectorized_linear.data(),
2938 sum_into_values);
2939}
2940
2941
2942
2943template <int n_components_, int dim, int spacedim, typename Number>
2944template <bool do_JxW, std::size_t stride_view>
2945inline void
2947 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2948 const EvaluationFlags::EvaluationFlags &integration_flags,
2949 const bool sum_into_values)
2950{
2951 // slow path with FEValues
2952 Assert(this->fe_values.get() != nullptr,
2953 ExcMessage(
2954 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2955 if (!sum_into_values)
2956 for (unsigned int i = 0; i < solution_values.size(); ++i)
2957 solution_values[i] = 0;
2958
2959 const std::size_t n_points = this->fe_values->get_quadrature().size();
2960
2961 if (integration_flags & EvaluationFlags::values)
2962 {
2963 AssertIndexRange(this->n_q_points, this->values.size() + 1);
2964 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2965 {
2966 for (unsigned int d = 0; d < n_components; ++d)
2967 if (this->nonzero_shape_function_component[i][d] &&
2968 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2969 for (unsigned int qb = 0, q = 0; q < n_points;
2970 ++qb, q += n_lanes_user_interface)
2971 for (unsigned int v = 0;
2972 v < n_lanes_user_interface && q + v < n_points;
2973 ++v)
2974 solution_values[i] +=
2975 this->fe_values->shape_value(i, q + v) *
2976 ETT::access(this->values[qb], v, d) *
2977 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2978 else if (this->nonzero_shape_function_component[i][d])
2979 for (unsigned int qb = 0, q = 0; q < n_points;
2980 ++qb, q += n_lanes_user_interface)
2981 for (unsigned int v = 0;
2982 v < n_lanes_user_interface && q + v < n_points;
2983 ++v)
2984 solution_values[i] +=
2985 this->fe_values->shape_value_component(i, q + v, d) *
2986 ETT::access(this->values[qb], v, d) *
2987 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2988 }
2989 }
2990
2991 if (integration_flags & EvaluationFlags::gradients)
2992 {
2993 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
2994 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2995 {
2996 for (unsigned int d = 0; d < n_components; ++d)
2997 if (this->nonzero_shape_function_component[i][d] &&
2998 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2999 for (unsigned int qb = 0, q = 0; q < n_points;
3000 ++qb, q += n_lanes_user_interface)
3001 for (unsigned int v = 0;
3002 v < n_lanes_user_interface && q + v < n_points;
3003 ++v)
3004 solution_values[i] +=
3005 this->fe_values->shape_grad(i, q + v) *
3006 ETT::access(this->gradients[qb], v, d) *
3007 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3008 else if (this->nonzero_shape_function_component[i][d])
3009 for (unsigned int qb = 0, q = 0; q < n_points;
3010 ++qb, q += n_lanes_user_interface)
3011 for (unsigned int v = 0;
3012 v < n_lanes_user_interface && q + v < n_points;
3013 ++v)
3014 solution_values[i] +=
3015 this->fe_values->shape_grad_component(i, q + v, d) *
3016 ETT::access(this->gradients[qb], v, d) *
3017 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3018 }
3019 }
3020}
3021
3022
3023
3024template <int n_components_, int dim, int spacedim, typename Number>
3025template <bool do_JxW, std::size_t stride_view>
3026void
3028 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3029 const EvaluationFlags::EvaluationFlags &integration_flags,
3030 const bool sum_into_values)
3031{
3032 if (this->must_reinitialize_pointers)
3033 internal_reinit_single_cell_state_mapping_info();
3034
3035 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3036
3037 if (this->n_q_points == 0 || // no evaluation points provided
3038 !((integration_flags & EvaluationFlags::values) ||
3039 (integration_flags &
3040 EvaluationFlags::gradients))) // no integration flags
3041 {
3042 if (!sum_into_values)
3043 for (unsigned int i = 0; i < solution_values.size(); ++i)
3044 solution_values[i] = 0;
3045 return;
3046 }
3047
3048 Assert(
3049 !do_JxW || this->JxW_ptr != nullptr,
3050 ExcMessage(
3051 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3052
3053 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3054 if (this->fast_path)
3055 {
3056 if (this->use_linear_path)
3057 integrate_fast<do_JxW, true>(solution_values,
3058 integration_flags,
3059 sum_into_values);
3060 else
3061 integrate_fast<do_JxW, false>(solution_values,
3062 integration_flags,
3063 sum_into_values);
3064 }
3065 else
3066 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3067}
3068
3069
3070
3071template <int n_components_, int dim, int spacedim, typename Number>
3074 const unsigned int point_index) const
3075{
3076 AssertIndexRange(point_index, this->n_q_points);
3077 Assert(this->normal_ptr != nullptr,
3079 ExcFEPointEvaluationAccessToUninitializedMappingField(
3080 "update_normal_vectors"));
3081 if (this->is_interior)
3082 return this->normal_ptr[point_index];
3083 else
3084 return -this->normal_ptr[point_index];
3085}
3086
3087
3088
3089template <int n_components_, int dim, int spacedim, typename Number>
3094 const bool is_interior,
3095 const unsigned int first_selected_component)
3096 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3097 mapping_info,
3098 fe,
3099 first_selected_component,
3100 is_interior)
3101{}
3102
3103
3104
3105template <int n_components_, int dim, int spacedim, typename Number>
3106inline void
3108 const unsigned int cell_index,
3109 const unsigned int face_number)
3110{
3111 this->current_cell_index = cell_index;
3112 this->current_face_number = face_number;
3113 this->must_reinitialize_pointers = false;
3114
3115 if (this->use_linear_path)
3116 this->template do_reinit<true, true>();
3117 else
3118 this->template do_reinit<true, false>();
3119}
3120
3121
3122
3123template <int n_components_, int dim, int spacedim, typename Number>
3124inline void
3126 const unsigned int face_index)
3127{
3128 this->current_cell_index = face_index;
3129 this->current_face_number =
3130 this->mapping_info->get_face_number(face_index, this->is_interior);
3131 this->must_reinitialize_pointers = false;
3132
3133 if (this->use_linear_path)
3134 this->template do_reinit<true, true>();
3135 else
3136 this->template do_reinit<true, false>();
3137}
3138
3139
3140
3141template <int n_components_, int dim, int spacedim, typename Number>
3142template <std::size_t stride_view>
3143void
3146 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3147{
3148 Assert(!this->must_reinitialize_pointers,
3149 ExcMessage("Object has not been reinitialized!"));
3150
3151 if (this->n_q_points == 0)
3152 return;
3153
3154 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3155
3156 if (!((evaluation_flags & EvaluationFlags::values) ||
3157 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3158 return;
3159
3160 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3161
3162 if (this->use_linear_path)
3163 do_evaluate<true>(solution_values, evaluation_flags);
3164 else
3165 do_evaluate<false>(solution_values, evaluation_flags);
3166}
3167
3168
3169
3170template <int n_components_, int dim, int spacedim, typename Number>
3171void
3173 const ArrayView<const ScalarNumber> &solution_values,
3174 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3175{
3176 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3177 solution_values.size()),
3178 evaluation_flags);
3179}
3180
3181
3182
3183template <int n_components_, int dim, int spacedim, typename Number>
3184template <bool is_linear, std::size_t stride_view>
3185void
3188 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3189{
3190 const unsigned int dofs_per_comp =
3191 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3192
3193 const ScalarNumber *input;
3194 if (stride_view == 1 && this->component_in_base_element == 0 &&
3195 (is_linear || this->renumber.empty()))
3196 input = solution_values.data();
3197 else
3198 {
3199 for (unsigned int comp = 0; comp < n_components; ++comp)
3200 {
3201 const std::size_t offset =
3202 (this->component_in_base_element + comp) * dofs_per_comp;
3203
3204 if (is_linear || this->renumber.empty())
3205 {
3206 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3207 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3208 solution_values[i + offset];
3209 }
3210 else
3211 {
3212 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3213 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3214 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3215 solution_values[renumber_ptr[i]];
3216 }
3217 }
3218 input = this->scratch_data_scalar.data();
3219 }
3220
3221 ScalarNumber *output =
3222 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3223
3224 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3225 template interpolate<true, false>(n_components,
3226 evaluation_flags,
3227 this->shape_info,
3228 input,
3229 output,
3230 this->current_face_number);
3231
3232 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3233}
3234
3235
3236
3237template <int n_components_, int dim, int spacedim, typename Number>
3238template <std::size_t stride_view>
3239void
3241 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3242 const EvaluationFlags::EvaluationFlags &integration_flags,
3243 const bool sum_into_values)
3244{
3245 Assert(!this->must_reinitialize_pointers,
3246 ExcMessage("Object has not been reinitialized!"));
3247
3248 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3249
3250 if (this->n_q_points == 0 || // no evaluation points provided
3251 !((integration_flags & EvaluationFlags::values) ||
3252 (integration_flags &
3253 EvaluationFlags::gradients))) // no integration flags
3254 {
3255 if (!sum_into_values)
3256 for (unsigned int i = 0; i < solution_values.size(); ++i)
3257 solution_values[i] = 0;
3258 return;
3259 }
3260
3261 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3262
3263 if (this->use_linear_path)
3264 do_integrate<true, true>(solution_values,
3265 integration_flags,
3266 sum_into_values);
3267 else
3268 do_integrate<true, false>(solution_values,
3269 integration_flags,
3270 sum_into_values);
3271}
3272
3273
3274
3275template <int n_components_, int dim, int spacedim, typename Number>
3276void
3278 const ArrayView<ScalarNumber> &solution_values,
3279 const EvaluationFlags::EvaluationFlags &integration_flags,
3280 const bool sum_into_values)
3281{
3282 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3283 solution_values.size()),
3284 integration_flags,
3285 sum_into_values);
3286}
3287
3288
3289
3290template <int n_components_, int dim, int spacedim, typename Number>
3291template <std::size_t stride_view>
3292void
3294 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3295 const EvaluationFlags::EvaluationFlags &integration_flags,
3296 const bool sum_into_values)
3297{
3298 Assert(!this->must_reinitialize_pointers,
3299 ExcMessage("Object has not been reinitialized!"));
3300
3301 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3302
3303 if (this->n_q_points == 0 || // no evaluation points provided
3304 !((integration_flags & EvaluationFlags::values) ||
3305 (integration_flags &
3306 EvaluationFlags::gradients))) // no integration flags
3307 {
3308 if (!sum_into_values)
3309 for (unsigned int i = 0; i < solution_values.size(); ++i)
3310 solution_values[i] = 0;
3311 return;
3312 }
3313
3314 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3315
3316 if (this->use_linear_path)
3317 do_integrate<false, true>(solution_values,
3318 integration_flags,
3319 sum_into_values);
3320 else
3321 do_integrate<false, false>(solution_values,
3322 integration_flags,
3323 sum_into_values);
3324}
3325
3326
3327
3328template <int n_components_, int dim, int spacedim, typename Number>
3329void
3331 const ArrayView<ScalarNumber> &solution_values,
3332 const EvaluationFlags::EvaluationFlags &integration_flags,
3333 const bool sum_into_values)
3334{
3335 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3336 solution_values.size()),
3337 integration_flags,
3338 sum_into_values);
3339}
3340
3341
3342
3343template <int n_components_, int dim, int spacedim, typename Number>
3344template <bool do_JxW, bool is_linear, std::size_t stride_view>
3345void
3347 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3348 const EvaluationFlags::EvaluationFlags &integration_flags,
3349 const bool sum_into_values)
3350{
3351 if (!sum_into_values && this->fe->n_components() > n_components)
3352 for (unsigned int i = 0; i < solution_values.size(); ++i)
3353 solution_values[i] = 0;
3354
3355 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3356 integration_flags,
3357 false);
3358
3359 ScalarNumber *input = this->scratch_data_scalar.begin();
3360
3361 if (stride_view == 1 && this->component_in_base_element == 0 &&
3362 (is_linear || this->renumber.empty()))
3363 {
3364 if (sum_into_values)
3365 internal::
3366 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3367 template interpolate<false, true>(n_components,
3368 integration_flags,
3369 this->shape_info,
3370 input,
3371 solution_values.data(),
3372 this->current_face_number);
3373 else
3374 internal::
3375 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3376 template interpolate<false, false>(n_components,
3377 integration_flags,
3378 this->shape_info,
3379 input,
3380 solution_values.data(),
3381 this->current_face_number);
3382 }
3383 else
3384 {
3385 const unsigned int dofs_per_comp_face =
3386 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3387
3388 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3389 ScalarNumber *output = input + size_input;
3390
3391 internal::
3392 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3393 template interpolate<false, false>(n_components,
3394 integration_flags,
3395 this->shape_info,
3396 input,
3397 output,
3398 this->current_face_number);
3399
3400 const unsigned int dofs_per_comp =
3401 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3402
3403 for (unsigned int comp = 0; comp < n_components; ++comp)
3404 {
3405 const std::size_t offset =
3406 (this->component_in_base_element + comp) * dofs_per_comp;
3407
3408 if (is_linear || this->renumber.empty())
3409 {
3410 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3411 if (sum_into_values)
3412 solution_values[i + offset] +=
3413 output[i + comp * dofs_per_comp];
3414 else
3415 solution_values[i + offset] =
3416 output[i + comp * dofs_per_comp];
3417 }
3418 else
3419 {
3420 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3421 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3422 if (sum_into_values)
3423 solution_values[renumber_ptr[i]] +=
3424 output[i + comp * dofs_per_comp];
3425 else
3426 solution_values[renumber_ptr[i]] =
3427 output[i + comp * dofs_per_comp];
3428 }
3429 }
3430 }
3431}
3432
3433
3434
3435template <int n_components_, int dim, int spacedim, typename Number>
3436template <int stride_face_dof>
3437void
3439 const ScalarNumber *face_dof_values,
3440 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3441{
3442 if (this->use_linear_path)
3443 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3444 evaluation_flags);
3445 else
3446 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3447 evaluation_flags);
3448}
3449
3450
3451
3452template <int n_components_, int dim, int spacedim, typename Number>
3453template <bool is_linear, int stride_face_dof>
3454inline void
3456 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3457 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3458{
3459 const scalar_value_type *face_dof_values_ptr;
3460 if constexpr (n_components == 1)
3461 face_dof_values_ptr = face_dof_values;
3462 else
3463 {
3464 const unsigned int dofs_per_comp_face =
3465 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3466 for (unsigned int comp = 0; comp < n_components; ++comp)
3467 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3468 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3469 stride_face_dof],
3470 comp,
3471 this->solution_renumbered[i]);
3472
3473 face_dof_values_ptr = this->solution_renumbered.data();
3474 }
3475
3476 constexpr int stride_face_dof_actual =
3477 n_components == 1 ? stride_face_dof : 1;
3478
3479 // loop over quadrature batches qb
3480 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3481
3482 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3483 {
3486
3487 if (evaluation_flags & EvaluationFlags::gradients)
3488 {
3489 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3490 is_linear ?
3492 dim - 1,
3495 2,
3496 stride_face_dof_actual>(face_dof_values_ptr,
3497 this->unit_point_faces_ptr[qb]) :
3499 dim - 1,
3502 2,
3503 false,
3504 stride_face_dof_actual>(this->shapes_faces.data() +
3505 qb * n_shapes,
3506 n_shapes,
3507 face_dof_values_ptr);
3508
3509 value = interpolated_value[dim - 1];
3510 // reorder derivative from tangential/normal derivatives into tensor
3511 // in physical coordinates
3512 if (this->current_face_number / 2 == 0)
3513 {
3514 gradient[0] = interpolated_value[dim];
3515 if (dim > 1)
3516 gradient[1] = interpolated_value[0];
3517 if (dim > 2)
3518 gradient[2] = interpolated_value[1];
3519 }
3520 else if (this->current_face_number / 2 == 1)
3521 {
3522 if (dim > 1)
3523 gradient[1] = interpolated_value[dim];
3524 if (dim == 3)
3525 {
3526 gradient[0] = interpolated_value[1];
3527 gradient[2] = interpolated_value[0];
3528 }
3529 else if (dim == 2)
3530 gradient[0] = interpolated_value[0];
3531 else
3533 }
3534 else if (this->current_face_number / 2 == 2)
3535 {
3536 if (dim > 2)
3537 {
3538 gradient[0] = interpolated_value[0];
3539 gradient[1] = interpolated_value[1];
3540 gradient[2] = interpolated_value[dim];
3541 }
3542 else
3544 }
3545 else
3547 }
3548 else
3549 {
3550 value = is_linear ?
3552 dim - 1,
3555 stride_face_dof_actual>(face_dof_values_ptr,
3556 this->unit_point_faces_ptr[qb]) :
3558 dim - 1,
3561 false,
3562 stride_face_dof_actual>(this->shapes_faces.data() +
3563 qb * n_shapes,
3564 n_shapes,
3565 face_dof_values_ptr);
3566 }
3567
3568 if (evaluation_flags & EvaluationFlags::values)
3569 {
3570 for (unsigned int v = 0, offset = qb * stride;
3571 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3572 ++v, ++offset)
3573 ETT::set_value(value, v, this->values[offset]);
3574 }
3575 if (evaluation_flags & EvaluationFlags::gradients)
3576 {
3577 Assert(this->update_flags & update_gradients ||
3578 this->update_flags & update_inverse_jacobians,
3580
3581 for (unsigned int v = 0, offset = qb * stride;
3582 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3583 ++v, ++offset)
3584 {
3585 unit_gradient_type unit_gradient;
3586 ETT::set_gradient(gradient, v, unit_gradient);
3587 this->gradients[offset] =
3588 this->cell_type <=
3591 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3593 this
3594 ->inverse_jacobian_ptr[this->cell_type <=
3596 GeometryType::affine ?
3597 0 :
3598 offset]
3599 .transpose(),
3600 unit_gradient);
3601 }
3602 }
3603 }
3604}
3605
3606
3607
3608template <int n_components_, int dim, int spacedim, typename Number>
3609template <int stride_face_dof>
3610void
3612 ScalarNumber *face_dof_values,
3613 const EvaluationFlags::EvaluationFlags &integration_flags,
3614 const bool sum_into_values)
3615{
3616 if (this->use_linear_path)
3617 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3618 integration_flags,
3619 sum_into_values);
3620 else
3621 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3622 integration_flags,
3623 sum_into_values);
3624}
3625
3626
3627
3628template <int n_components_, int dim, int spacedim, typename Number>
3629template <bool do_JxW, bool is_linear, int stride_face_dof>
3630inline void
3633 ScalarNumber *face_dof_values,
3634 const EvaluationFlags::EvaluationFlags &integration_flags,
3635 const bool sum_into_values)
3636{
3637 // zero out lanes of incomplete last quadrature point batch
3638 if constexpr (stride == 1)
3639 if (const unsigned int n_filled_lanes =
3640 this->n_q_points_scalar & (n_lanes_internal - 1);
3641 n_filled_lanes > 0)
3642 {
3643 if (integration_flags & EvaluationFlags::values)
3644 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3645 ETT::set_zero_value(this->values.back(), v);
3646 if (integration_flags & EvaluationFlags::gradients)
3647 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3648 ETT::set_zero_gradient(this->gradients.back(), v);
3649 }
3650
3651 std::array<vectorized_value_type,
3652 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3653 solution_values_vectorized_linear = {};
3654
3655 // loop over quadrature batches qb
3656 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3657
3658 const bool cartesian_cell =
3660 const bool affine_cell =
3662 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3663 {
3664 vectorized_value_type value = {};
3666
3667 if (integration_flags & EvaluationFlags::values)
3668 for (unsigned int v = 0, offset = qb * stride;
3669 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3670 ++v, ++offset)
3671 ETT::get_value(value,
3672 v,
3673 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3674 this->values[offset]);
3675
3676 if (integration_flags & EvaluationFlags::gradients)
3677 for (unsigned int v = 0, offset = qb * stride;
3678 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3679 ++v, ++offset)
3680 {
3681 const auto grad_w =
3682 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3683 this->gradients[offset];
3684 ETT::get_gradient(
3685 gradient,
3686 v,
3687 cartesian_cell ?
3688 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3689 grad_w) :
3691 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3692 grad_w));
3693 }
3694
3695 if (integration_flags & EvaluationFlags::gradients)
3696 {
3697 std::array<vectorized_value_type, 2> value_face = {};
3698 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3699
3700 value_face[0] = value;
3701 // fill derivative in physical coordinates into tangential/normal
3702 // derivatives
3703 if (this->current_face_number / 2 == 0)
3704 {
3705 value_face[1] = gradient[0];
3706 if (dim > 1)
3707 gradient_in_face[0] = gradient[1];
3708 if (dim > 2)
3709 gradient_in_face[1] = gradient[2];
3710 }
3711 else if (this->current_face_number / 2 == 1)
3712 {
3713 if (dim > 1)
3714 value_face[1] = gradient[1];
3715 if (dim == 3)
3716 {
3717 gradient_in_face[0] = gradient[2];
3718 gradient_in_face[1] = gradient[0];
3719 }
3720 else if (dim == 2)
3721 gradient_in_face[0] = gradient[0];
3722 else
3724 }
3725 else if (this->current_face_number / 2 == 2)
3726 {
3727 if (dim > 2)
3728 {
3729 value_face[1] = gradient[2];
3730 gradient_in_face[0] = gradient[0];
3731 gradient_in_face[1] = gradient[1];
3732 }
3733 else
3735 }
3736 else
3738
3740 is_linear,
3741 dim - 1,
3744 2>(this->shapes_faces.data() + qb * n_shapes,
3745 n_shapes,
3746 value_face.data(),
3747 gradient_in_face,
3748 is_linear ? solution_values_vectorized_linear.data() :
3749 this->solution_renumbered_vectorized.data(),
3750 this->unit_point_faces_ptr[qb],
3751 qb != 0);
3752 }
3753 else
3755 dim - 1,
3758 this->shapes_faces.data() + qb * n_shapes,
3759 n_shapes,
3760 value,
3761 is_linear ? solution_values_vectorized_linear.data() :
3762 this->solution_renumbered_vectorized.data(),
3763 this->unit_point_faces_ptr[qb],
3764 qb != 0);
3765 }
3766
3767 const unsigned int dofs_per_comp_face =
3768 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3769
3770 for (unsigned int comp = 0; comp < n_components; ++comp)
3771 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3772 if (sum_into_values)
3773 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3774 stride_face_dof] +=
3775 ETT::sum_value(comp,
3776 is_linear ?
3777 *(solution_values_vectorized_linear.data() + i) :
3778 this->solution_renumbered_vectorized[i]);
3779 else
3780 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3781 ETT::sum_value(comp,
3782 is_linear ?
3783 *(solution_values_vectorized_linear.data() + i) :
3784 this->solution_renumbered_vectorized[i]);
3785}
3786
3787
3788
3789template <int n_components_, int dim, int spacedim, typename Number>
3792 const unsigned int point_index) const
3793{
3794 AssertIndexRange(point_index, this->n_q_points);
3795 Assert(this->normal_ptr != nullptr,
3797 ExcFEPointEvaluationAccessToUninitializedMappingField(
3798 "update_normal_vectors"));
3799 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3800 {
3802 for (unsigned int d = 0; d < dim; ++d)
3803 normal[d] =
3804 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3805 0);
3806 if (this->is_interior)
3807 return normal;
3808 else
3809 return -normal;
3810 }
3811 else
3812 {
3813 if (this->is_interior)
3814 return this->normal_ptr[point_index];
3815 else
3816 return -(this->normal_ptr[point_index]);
3817 }
3818}
3819
3821
3822#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 internal_reinit_single_cell_state_mapping_info()
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_unit_gradient_type &gradient)
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr std::size_t n_lanes_user_interface
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
Abstract base class for mapping classes.
Definition mapping.h: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)