Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07: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
707
708public:
716 const value_type &
717 get_value(const unsigned int point_index) const;
718
727 void
728 submit_value(const value_type &value, const unsigned int point_index);
729
739 const gradient_type &
740 get_gradient(const unsigned int point_index) const;
741
750 void
751 submit_gradient(const gradient_type &, const unsigned int point_index);
752
759 jacobian(const unsigned int point_index) const;
760
768 inverse_jacobian(const unsigned int point_index) const;
769
775 Number
776 JxW(const unsigned int point_index) const;
777
784 DEAL_II_DEPRECATED_EARLY Point<spacedim, Number>
785 real_point(const unsigned int point_index) const;
786
792 quadrature_point(const unsigned int point_index) const;
793
799 unit_point(const unsigned int point_index) const;
800
810
818
822 unsigned int
824
825protected:
826 static constexpr std::size_t n_lanes_user_interface =
828 static constexpr std::size_t n_lanes_internal =
830 static constexpr std::size_t stride =
832
841 void
842 setup(const unsigned int first_selected_component);
843
849 template <bool is_face, bool is_linear>
850 void
852
856 const unsigned int n_q_batches;
857
861 const unsigned int n_q_points;
862
866 const unsigned int n_q_points_scalar;
867
872
877
882 std::vector<Polynomials::Polynomial<double>> poly;
883
888
893 std::vector<unsigned int> renumber;
894
902 std::vector<scalar_value_type> solution_renumbered;
903
911
916
920 std::vector<value_type> values;
921
925 std::vector<gradient_type> gradients;
926
932
938
944
950
956
962
968 const Number *JxW_ptr;
969
974
979 unsigned int dofs_per_component;
980
987
992
997
1003 std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
1004
1009
1013 std::shared_ptr<FEValues<dim, spacedim>> fe_values;
1014
1018 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1020
1026
1031
1036
1041
1046 boost::signals2::connection connection_is_reinitialized;
1047
1052
1058
1064
1065 const bool is_interior;
1066};
1067
1068
1069
1100template <int n_components_,
1101 int dim,
1102 int spacedim = dim,
1103 typename Number = double>
1105 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1106{
1107public:
1108 static constexpr unsigned int dimension = dim;
1109 static constexpr unsigned int n_components = n_components_;
1110
1111 using number_type = Number;
1112
1115 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1117 using ETT = typename internal::FEPointEvaluation::
1118 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1119 using value_type = typename ETT::value_type;
1120 using scalar_value_type = typename ETT::scalar_value_type;
1121 using vectorized_value_type = typename ETT::vectorized_value_type;
1122 using unit_gradient_type = typename ETT::unit_gradient_type;
1123 using gradient_type = typename ETT::real_gradient_type;
1125 typename ETT::interface_vectorized_unit_gradient_type;
1126
1148 const unsigned int first_selected_component = 0);
1149
1169 const unsigned int first_selected_component = 0);
1170
1182 void
1184 const ArrayView<const Point<dim>> &unit_points);
1185
1190 void
1191 reinit();
1192
1197 void
1198 reinit(const unsigned int cell_index);
1199
1200
1212 template <std::size_t stride_view>
1213 void
1214 evaluate(
1216 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1217
1229 void
1230 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1231 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1232
1255 template <std::size_t stride_view>
1256 void
1258 const EvaluationFlags::EvaluationFlags &integration_flags,
1259 const bool sum_into_values = false);
1260
1283 void
1284 integrate(const ArrayView<ScalarNumber> &solution_values,
1285 const EvaluationFlags::EvaluationFlags &integration_flags,
1286 const bool sum_into_values = false);
1287
1314 template <std::size_t stride_view>
1315 void
1317 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1318 const EvaluationFlags::EvaluationFlags &integration_flags,
1319 const bool sum_into_values = false);
1320
1347 void
1348 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1349 const EvaluationFlags::EvaluationFlags &integration_flags,
1350 const bool sum_into_values = false);
1351
1358 normal_vector(const unsigned int point_index) const;
1359
1360private:
1361 static constexpr std::size_t n_lanes_user_interface =
1363 static constexpr std::size_t n_lanes_internal =
1365 static constexpr std::size_t stride =
1367
1372 template <bool is_linear, std::size_t stride_view>
1373 void
1376
1381 template <bool is_linear, std::size_t stride_view>
1382 void
1385 const EvaluationFlags::EvaluationFlags &evaluation_flags,
1386 const unsigned int n_shapes,
1387 const unsigned int qb,
1388 vectorized_value_type &value,
1390
1394 template <bool is_linear, std::size_t stride_view>
1395 void
1398 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1399
1403 template <std::size_t stride_view>
1404 void
1407 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1408
1414 template <bool is_linear>
1415 void
1417 const EvaluationFlags::EvaluationFlags &integration_flags,
1418 const unsigned int n_shapes,
1419 const unsigned int qb,
1420 const vectorized_value_type value,
1422 vectorized_value_type *solution_values_vectorized_linear);
1423
1429 template <bool is_linear, std::size_t stride_view>
1430 void
1432 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1433 vectorized_value_type *solution_values_vectorized_linear,
1434 const bool sum_into_values);
1435
1439 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1440 void
1442 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1443 const EvaluationFlags::EvaluationFlags &integration_flags,
1444 const bool sum_into_values);
1445
1449 template <bool do_JxW, std::size_t stride_view>
1450 void
1452 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1453 const EvaluationFlags::EvaluationFlags &integration_flags,
1454 const bool sum_into_values);
1455
1459 template <bool do_JxW, std::size_t stride_view>
1460 void
1462 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1463 const EvaluationFlags::EvaluationFlags &integration_flags,
1464 const bool sum_into_values);
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 , is_reinitialized(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 , is_reinitialized(false)
1803 , is_interior(is_interior)
1804{
1805 setup(first_selected_component);
1806 connection_is_reinitialized = mapping_info.connect_is_reinitialized(
1807 [this]() { this->is_reinitialized = false; });
1808}
1809
1810
1811
1812template <int n_components_, int dim, int spacedim, typename Number>
1816 : n_q_batches(other.n_q_batches)
1817 , n_q_points(other.n_q_points)
1818 , n_q_points_scalar(other.n_q_points_scalar)
1819 , mapping(other.mapping)
1820 , fe(other.fe)
1821 , poly(other.poly)
1822 , use_linear_path(other.use_linear_path)
1823 , renumber(other.renumber)
1824 , solution_renumbered(other.solution_renumbered)
1825 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1826 , values(other.values)
1827 , gradients(other.gradients)
1828 , dofs_per_component(other.dofs_per_component)
1829 , dofs_per_component_face(other.dofs_per_component_face)
1830 , component_in_base_element(other.component_in_base_element)
1831 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1832 , update_flags(other.update_flags)
1833 , fe_values(other.fe_values)
1834 , mapping_info_on_the_fly(
1835 other.mapping_info_on_the_fly ?
1837 *mapping,
1838 update_flags) :
1839 nullptr)
1840 , mapping_info(other.mapping_info)
1841 , current_cell_index(other.current_cell_index)
1842 , current_face_number(other.current_face_number)
1843 , fast_path(other.fast_path)
1844 , is_reinitialized(false)
1845 , is_interior(other.is_interior)
1846{
1847 connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1848 [this]() { this->is_reinitialized = false; });
1849}
1850
1851
1852
1853template <int n_components_, int dim, int spacedim, typename Number>
1856 : n_q_batches(other.n_q_batches)
1857 , n_q_points(other.n_q_points)
1858 , n_q_points_scalar(other.n_q_points_scalar)
1859 , mapping(other.mapping)
1860 , fe(other.fe)
1861 , poly(other.poly)
1862 , use_linear_path(other.use_linear_path)
1863 , renumber(other.renumber)
1864 , solution_renumbered(other.solution_renumbered)
1865 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1866 , values(other.values)
1867 , gradients(other.gradients)
1868 , dofs_per_component(other.dofs_per_component)
1869 , dofs_per_component_face(other.dofs_per_component_face)
1870 , component_in_base_element(other.component_in_base_element)
1871 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1872 , update_flags(other.update_flags)
1873 , fe_values(other.fe_values)
1874 , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1875 , mapping_info(other.mapping_info)
1876 , current_cell_index(other.current_cell_index)
1877 , current_face_number(other.current_face_number)
1878 , fast_path(other.fast_path)
1879 , is_reinitialized(false)
1880 , is_interior(other.is_interior)
1881{
1882 connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1883 [this]() { this->is_reinitialized = false; });
1884}
1885
1886
1887
1888template <int n_components_, int dim, int spacedim, typename Number>
1891{
1892 connection_is_reinitialized.disconnect();
1893}
1894
1895
1896
1897template <int n_components_, int dim, int spacedim, typename Number>
1898void
1900 const unsigned int first_selected_component)
1901{
1902 AssertIndexRange(first_selected_component + n_components,
1903 fe->n_components() + 1);
1904
1905 shapes.reserve(100);
1906
1907 bool same_base_element = true;
1908 unsigned int base_element_number = 0;
1909 component_in_base_element = 0;
1910 unsigned int component = 0;
1911 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1912 if (component + fe->element_multiplicity(base_element_number) >
1913 first_selected_component)
1914 {
1915 if (first_selected_component + n_components >
1916 component + fe->element_multiplicity(base_element_number))
1917 same_base_element = false;
1918 component_in_base_element = first_selected_component - component;
1919 break;
1920 }
1921 else
1922 component += fe->element_multiplicity(base_element_number);
1923
1926 *fe, base_element_number) &&
1927 same_base_element)
1928 {
1929 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1930 renumber = shape_info.lexicographic_numbering;
1931 dofs_per_component = shape_info.dofs_per_component_on_cell;
1932 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1934 fe->base_element(base_element_number));
1935
1936 bool is_lexicographic = true;
1937 for (unsigned int i = 0; i < renumber.size(); ++i)
1938 if (i != renumber[i])
1939 is_lexicographic = false;
1940
1941 if (is_lexicographic)
1942 renumber.clear();
1943
1944 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1945 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1946 poly[1].value(1.) == 1.) &&
1947 (fe->n_components() == n_components);
1948
1949 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1950 const unsigned int size_cell = dofs_per_component * n_components;
1951 scratch_data_scalar.resize(size_face + size_cell);
1952
1953 solution_renumbered.resize(dofs_per_component);
1954 solution_renumbered_vectorized.resize(dofs_per_component);
1955
1956 fast_path = true;
1957 }
1958 else
1959 {
1960 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1961 for (unsigned int d = 0; d < n_components; ++d)
1962 {
1963 const unsigned int component = first_selected_component + d;
1964 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1965 {
1966 const bool is_primitive =
1967 fe->is_primitive() || fe->is_primitive(i);
1968 if (is_primitive)
1969 nonzero_shape_function_component[i][d] =
1970 (component == fe->system_to_component_index(i).first);
1971 else
1972 nonzero_shape_function_component[i][d] =
1973 (fe->get_nonzero_components(i)[component] == true);
1974 }
1975 }
1976
1977 fast_path = false;
1978 }
1979}
1980
1981
1982
1983template <int n_components_, int dim, int spacedim, typename Number>
1984template <bool is_face, bool is_linear>
1985inline void
1987{
1988 const unsigned int geometry_index =
1989 mapping_info->template compute_geometry_index_offset<is_face>(
1990 current_cell_index, current_face_number);
1991
1992 cell_type = mapping_info->get_cell_type(geometry_index);
1993
1994 const_cast<unsigned int &>(n_q_points_scalar) =
1995 mapping_info->get_n_q_points_unvectorized(geometry_index);
1996
1997 // round up n_q_points_scalar / n_lanes_internal
1998 const_cast<unsigned int &>(n_q_batches) =
1999 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
2000
2001 const unsigned int n_q_points_before = n_q_points;
2002
2003 const_cast<unsigned int &>(n_q_points) =
2004 (stride == 1) ? n_q_batches : n_q_points_scalar;
2005
2006 if (n_q_points != n_q_points_before)
2007 {
2008 if (update_flags & update_values)
2009 values.resize(n_q_points);
2010 if (update_flags & update_gradients)
2011 gradients.resize(n_q_points);
2012 }
2013
2014 if (n_q_points == 0)
2015 {
2016 is_reinitialized = true;
2017 return;
2018 }
2019
2020 // set unit point pointer
2021 const unsigned int unit_point_offset =
2022 mapping_info->compute_unit_point_index_offset(geometry_index);
2023
2024 if (is_face)
2025 unit_point_faces_ptr =
2026 mapping_info->get_unit_point_faces(unit_point_offset);
2027 else
2028 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2029
2030 // set data pointers
2031 const unsigned int data_offset =
2032 mapping_info->compute_data_index_offset(geometry_index);
2033 const unsigned int compressed_data_offset =
2034 mapping_info->compute_compressed_data_index_offset(geometry_index);
2035#ifdef DEBUG
2036 const UpdateFlags update_flags_mapping =
2037 mapping_info->get_update_flags_mapping();
2038 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2039 real_point_ptr = mapping_info->get_real_point(data_offset);
2040 if (update_flags_mapping & UpdateFlags::update_jacobians)
2041 jacobian_ptr =
2042 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2043 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2044 inverse_jacobian_ptr =
2045 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2046 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2047 normal_ptr = mapping_info->get_normal_vector(data_offset);
2048 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2049 JxW_ptr = mapping_info->get_JxW(data_offset);
2050#else
2051 real_point_ptr = mapping_info->get_real_point(data_offset);
2052 jacobian_ptr =
2053 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2054 inverse_jacobian_ptr =
2055 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2056 normal_ptr = mapping_info->get_normal_vector(data_offset);
2057 JxW_ptr = mapping_info->get_JxW(data_offset);
2058#endif
2059
2060 if (!is_linear && fast_path)
2061 {
2062 const std::size_t n_shapes = poly.size();
2063 if (is_face)
2064 shapes_faces.resize_fast(n_q_batches * n_shapes);
2065 else
2066 shapes.resize_fast(n_q_batches * n_shapes);
2067
2068 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2069 if (is_face)
2070 {
2071 if (dim > 1)
2072 {
2074 shapes_faces.data() + qb * n_shapes,
2075 poly,
2076 unit_point_faces_ptr[qb],
2077 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2078 }
2079 }
2080 else
2081 {
2082 if (update_flags & UpdateFlags::update_gradients)
2083 {
2084 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2085 poly,
2086 unit_point_ptr[qb],
2087 1);
2088 }
2089 else if (qb + 1 < n_q_batches)
2090 {
2091 // Use function with reduced overhead to compute for two
2092 // points at once
2094 shapes.data() + qb * n_shapes,
2095 poly,
2096 unit_point_ptr[qb],
2097 unit_point_ptr[qb + 1]);
2098 ++qb;
2099 }
2100 else
2101 {
2102 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2103 poly,
2104 unit_point_ptr[qb],
2105 0);
2106 }
2107 }
2108 }
2109
2110 is_reinitialized = true;
2111}
2112
2113
2114
2115template <int n_components_, int dim, int spacedim, typename Number>
2116inline const typename FEPointEvaluationBase<n_components_,
2117 dim,
2118 spacedim,
2119 Number>::value_type &
2121 const unsigned int point_index) const
2122{
2123 AssertIndexRange(point_index, values.size());
2124 return values[point_index];
2125}
2126
2127
2128
2129template <int n_components_, int dim, int spacedim, typename Number>
2130inline const typename FEPointEvaluationBase<n_components_,
2131 dim,
2132 spacedim,
2133 Number>::gradient_type &
2135 const unsigned int point_index) const
2136{
2137 AssertIndexRange(point_index, gradients.size());
2138 return gradients[point_index];
2139}
2140
2141
2142
2143template <int n_components_, int dim, int spacedim, typename Number>
2144inline void
2146 const value_type &value,
2147 const unsigned int point_index)
2148{
2149 AssertIndexRange(point_index, n_q_points);
2150 values[point_index] = value;
2151}
2152
2153
2154
2155template <int n_components_, int dim, int spacedim, typename Number>
2156inline void
2158 const gradient_type &gradient,
2159 const unsigned int point_index)
2160{
2161 AssertIndexRange(point_index, n_q_points);
2162 gradients[point_index] = gradient;
2163}
2164
2165
2166
2167template <int n_components_, int dim, int spacedim, typename Number>
2170 const unsigned int point_index) const
2171{
2172 AssertIndexRange(point_index, n_q_points);
2173 Assert(jacobian_ptr != nullptr,
2175 ExcFEPointEvaluationAccessToUninitializedMappingField(
2176 "update_jacobians"));
2177 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2178 GeometryType::affine ?
2179 0 :
2180 point_index];
2181}
2182
2183
2184
2185template <int n_components_, int dim, int spacedim, typename Number>
2188 const unsigned int point_index) const
2189{
2190 AssertIndexRange(point_index, n_q_points);
2191 Assert(inverse_jacobian_ptr != nullptr,
2193 ExcFEPointEvaluationAccessToUninitializedMappingField(
2194 "update_inverse_jacobians"));
2195 return inverse_jacobian_ptr
2196 [cell_type <=
2198 0 :
2199 point_index];
2200}
2201
2202
2203
2204template <int n_components_, int dim, int spacedim, typename Number>
2205inline Number
2207 const unsigned int point_index) const
2208{
2209 AssertIndexRange(point_index, n_q_points);
2210 Assert(JxW_ptr != nullptr,
2212 ExcFEPointEvaluationAccessToUninitializedMappingField(
2213 "update_JxW_values"));
2214 return JxW_ptr[point_index];
2215}
2216
2217
2218
2219template <int n_components_, int dim, int spacedim, typename Number>
2222 const unsigned int point_index) const
2223{
2224 return quadrature_point(point_index);
2225 ;
2226}
2227
2228
2229
2230template <int n_components_, int dim, int spacedim, typename Number>
2233 const unsigned int point_index) const
2234{
2235 AssertIndexRange(point_index, n_q_points);
2236 Assert(real_point_ptr != nullptr,
2238 ExcFEPointEvaluationAccessToUninitializedMappingField(
2239 "update_quadrature_points"));
2240 return real_point_ptr[point_index];
2241}
2242
2243
2244
2245template <int n_components_, int dim, int spacedim, typename Number>
2246inline Point<dim, Number>
2248 const unsigned int point_index) const
2249{
2250 AssertIndexRange(point_index, n_q_points);
2251 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2252 Point<dim, Number> unit_point;
2253 for (unsigned int d = 0; d < dim; ++d)
2255 unit_point_ptr[point_index / stride][d], point_index % stride);
2256 return unit_point;
2257}
2258
2259
2260
2261template <int n_components_, int dim, int spacedim, typename Number>
2268
2269
2270
2271template <int n_components_, int dim, int spacedim, typename Number>
2275 const unsigned int first_selected_component)
2276 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2277 mapping_info,
2278 fe,
2279 first_selected_component)
2280{}
2281
2282
2283
2284template <int n_components_, int dim, int spacedim, typename Number>
2286 const Mapping<dim, spacedim> &mapping,
2288 const UpdateFlags update_flags,
2289 const unsigned int first_selected_component)
2290 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2291 mapping,
2292 fe,
2293 update_flags,
2294 first_selected_component)
2295{}
2296
2297
2298
2299template <int n_components_, int dim, int spacedim, typename Number>
2300inline void
2302{
2303 this->current_cell_index = numbers::invalid_unsigned_int;
2304 this->current_face_number = numbers::invalid_unsigned_int;
2305
2306 if (this->use_linear_path)
2307 this->template do_reinit<false, true>();
2308 else
2309 this->template do_reinit<false, false>();
2310}
2311
2312
2313
2314template <int n_components_, int dim, int spacedim, typename Number>
2315inline void
2318 const ArrayView<const Point<dim>> &unit_points)
2319{
2320 // reinit is only allowed for mapping computation on the fly
2321 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2323
2324 this->mapping_info->reinit(cell, unit_points);
2325
2326 if (!this->fast_path)
2327 {
2328 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2329 *this->mapping,
2330 *this->fe,
2332 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2333 this->update_flags);
2334 this->fe_values->reinit(cell);
2335 }
2336
2337 if (this->use_linear_path)
2338 this->template do_reinit<false, true>();
2339 else
2340 this->template do_reinit<false, false>();
2341}
2342
2343
2344
2345template <int n_components_, int dim, int spacedim, typename Number>
2346inline void
2348 const unsigned int cell_index)
2349{
2350 this->current_cell_index = cell_index;
2351 this->current_face_number = numbers::invalid_unsigned_int;
2352
2353 if (this->use_linear_path)
2354 this->template do_reinit<false, true>();
2355 else
2356 this->template do_reinit<false, false>();
2357
2358 if (!this->fast_path)
2359 {
2360 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2361
2362 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2363 for (unsigned int d = 0; d < dim; ++d)
2364 unit_points[v][d] =
2365 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2366
2367 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2368 *this->mapping,
2369 *this->fe,
2371 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2372 this->update_flags);
2373
2374 this->fe_values->reinit(
2375 this->mapping_info->get_cell_iterator(this->current_cell_index));
2376 }
2377}
2378
2379
2380
2381template <int n_components_, int dim, int spacedim, typename Number>
2382template <std::size_t stride_view>
2383void
2386 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2387{
2388 if (!this->is_reinitialized)
2389 reinit();
2390
2391 if (this->n_q_points == 0)
2392 return;
2393
2394 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2395
2396 if (!((evaluation_flags & EvaluationFlags::values) ||
2397 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2398 return;
2399
2400 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2401 if (this->fast_path)
2402 {
2403 if (this->use_linear_path)
2404 evaluate_fast<true>(solution_values, evaluation_flags);
2405 else
2406 evaluate_fast<false>(solution_values, evaluation_flags);
2407 }
2408 else
2409 evaluate_slow(solution_values, evaluation_flags);
2410}
2411
2412
2413
2414template <int n_components_, int dim, int spacedim, typename Number>
2415void
2417 const ArrayView<const ScalarNumber> &solution_values,
2418 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2419{
2420 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2421 solution_values.size()),
2422 evaluation_flags);
2423}
2424
2425
2426
2427template <int n_components_, int dim, int spacedim, typename Number>
2428template <std::size_t stride_view>
2429void
2431 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2432 const EvaluationFlags::EvaluationFlags &integration_flags,
2433 const bool sum_into_values)
2434{
2435 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2436}
2437
2438
2439
2440template <int n_components_, int dim, int spacedim, typename Number>
2441void
2443 const ArrayView<ScalarNumber> &solution_values,
2444 const EvaluationFlags::EvaluationFlags &integration_flags,
2445 const bool sum_into_values)
2446{
2447 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2448 solution_values.size()),
2449 integration_flags,
2450 sum_into_values);
2451}
2452
2453
2454
2455template <int n_components_, int dim, int spacedim, typename Number>
2457 scalar_value_type
2459 const
2460{
2461 value_type return_value = {};
2462
2463 for (const auto point_index : this->quadrature_point_indices())
2464 return_value += values[point_index] * this->JxW(point_index);
2465
2466 return ETT::sum_value(return_value);
2467}
2468
2469
2470
2471template <int n_components_, int dim, int spacedim, typename Number>
2472unsigned int
2475{
2476 Assert(stride == 1,
2477 ExcMessage(
2478 "Calling this function only makes sense in fully vectorized mode."));
2479 if (q == n_q_batches - 1)
2480 {
2481 const unsigned int n_filled_lanes =
2482 n_q_points_scalar & (n_lanes_user_interface - 1);
2483 if (n_filled_lanes == 0)
2484 return n_lanes_user_interface;
2485 else
2486 return n_filled_lanes;
2487 }
2488 else
2489 return n_lanes_user_interface;
2490}
2491
2492
2493
2494template <int n_components_, int dim, int spacedim, typename Number>
2495template <std::size_t stride_view>
2496void
2498 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2499 const EvaluationFlags::EvaluationFlags &integration_flags,
2500 const bool sum_into_values)
2501{
2502 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2503}
2504
2505
2506
2507template <int n_components_, int dim, int spacedim, typename Number>
2508void
2510 const ArrayView<ScalarNumber> &solution_values,
2511 const EvaluationFlags::EvaluationFlags &integration_flags,
2512 const bool sum_into_values)
2513{
2514 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2515 solution_values.size()),
2516 integration_flags,
2517 sum_into_values);
2518}
2519
2520
2521
2522template <int n_components_, int dim, int spacedim, typename Number>
2523template <bool is_linear, std::size_t stride_view>
2524inline void
2527{
2528 const unsigned int dofs_per_comp =
2529 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2530
2531 for (unsigned int comp = 0; comp < n_components; ++comp)
2532 {
2533 const std::size_t offset =
2534 (this->component_in_base_element + comp) * dofs_per_comp;
2535
2536 if ((is_linear && n_components == 1) || this->renumber.empty())
2537 {
2538 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2539 ETT::read_value(solution_values[i + offset],
2540 comp,
2541 this->solution_renumbered[i]);
2542 }
2543 else
2544 {
2545 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2546 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2547 ETT::read_value(solution_values[renumber_ptr[i]],
2548 comp,
2549 this->solution_renumbered[i]);
2550 }
2551 }
2552}
2553
2554
2555
2556template <int n_components_, int dim, int spacedim, typename Number>
2557template <bool is_linear, std::size_t stride_view>
2558inline void
2561 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2562 const unsigned int n_shapes,
2563 const unsigned int qb,
2564 vectorized_value_type &value,
2566{
2567 if (evaluation_flags & EvaluationFlags::gradients)
2568 {
2569 std::array<vectorized_value_type, dim + 1> result;
2570 if constexpr (is_linear)
2571 {
2572 if constexpr (n_components == 1)
2573 result =
2575 dim,
2578 1,
2579 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2580 else
2581 result =
2583 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2584 }
2585 else
2587 dim,
2590 1,
2591 false>(this->shapes.data() + qb * n_shapes,
2592 n_shapes,
2593 this->solution_renumbered.data());
2594 gradient[0] = result[0];
2595 if (dim > 1)
2596 gradient[1] = result[1];
2597 if (dim > 2)
2598 gradient[2] = result[2];
2599 value = result[dim];
2600 }
2601 else
2602 {
2603 if constexpr (is_linear)
2604 {
2605 if constexpr (n_components == 1)
2607 dim,
2610 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2611 else
2613 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2614 }
2615 else
2616 value =
2620 false>(
2621 this->shapes.data() + qb * n_shapes,
2622 n_shapes,
2623 this->solution_renumbered.data());
2624 }
2625}
2626
2627
2628
2629template <int n_components_, int dim, int spacedim, typename Number>
2630template <bool is_linear, std::size_t stride_view>
2631inline void
2634 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2635{
2636 if (!(is_linear && n_components == 1))
2637 prepare_evaluate_fast<is_linear>(solution_values);
2638
2639 // loop over quadrature batches qb
2640 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2641
2642 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2643 {
2646
2647 compute_evaluate_fast<is_linear>(
2648 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2649
2650 if (evaluation_flags & EvaluationFlags::values)
2651 {
2652 for (unsigned int v = 0, offset = qb * stride;
2653 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2654 ++v, ++offset)
2655 ETT::set_value(value, v, this->values[offset]);
2656 }
2657 if (evaluation_flags & EvaluationFlags::gradients)
2658 {
2659 Assert(this->update_flags & update_gradients ||
2660 this->update_flags & update_inverse_jacobians,
2662
2663 for (unsigned int v = 0, offset = qb * stride;
2664 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2665 ++v, ++offset)
2666 {
2667 unit_gradient_type unit_gradient;
2668 ETT::set_gradient(gradient, v, unit_gradient);
2669 this->gradients[offset] =
2670 this->cell_type <=
2673 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2675 this
2676 ->inverse_jacobian_ptr[this->cell_type <=
2678 GeometryType::affine ?
2679 0 :
2680 offset]
2681 .transpose(),
2682 unit_gradient);
2683 }
2684 }
2685 }
2686}
2687
2688
2689
2690template <int n_components_, int dim, int spacedim, typename Number>
2691template <std::size_t stride_view>
2692inline void
2695 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2696{
2697 // slow path with FEValues
2698 Assert(this->fe_values.get() != nullptr,
2699 ExcMessage(
2700 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2701
2702 const std::size_t n_points = this->fe_values->get_quadrature().size();
2703
2704 if (evaluation_flags & EvaluationFlags::values)
2705 {
2706 this->values.resize(this->n_q_points);
2707 std::fill(this->values.begin(), this->values.end(), value_type());
2708 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2709 {
2710 const ScalarNumber value = solution_values[i];
2711 for (unsigned int d = 0; d < n_components; ++d)
2712 if (this->nonzero_shape_function_component[i][d] &&
2713 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2714 for (unsigned int qb = 0, q = 0; q < n_points;
2715 ++qb, q += n_lanes_user_interface)
2716 for (unsigned int v = 0;
2717 v < n_lanes_user_interface && q + v < n_points;
2718 ++v)
2719 ETT::access(this->values[qb],
2720 v,
2721 d,
2722 this->fe_values->shape_value(i, q + v) * value);
2723 else if (this->nonzero_shape_function_component[i][d])
2724 for (unsigned int qb = 0, q = 0; q < n_points;
2725 ++qb, q += n_lanes_user_interface)
2726 for (unsigned int v = 0;
2727 v < n_lanes_user_interface && q + v < n_points;
2728 ++v)
2729 ETT::access(this->values[qb],
2730 v,
2731 d,
2732 this->fe_values->shape_value_component(i,
2733 q + v,
2734 d) *
2735 value);
2736 }
2737 }
2738
2739 if (evaluation_flags & EvaluationFlags::gradients)
2740 {
2741 this->gradients.resize(this->n_q_points);
2742 std::fill(this->gradients.begin(),
2743 this->gradients.end(),
2744 gradient_type());
2745 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2746 {
2747 const ScalarNumber value = solution_values[i];
2748 for (unsigned int d = 0; d < n_components; ++d)
2749 if (this->nonzero_shape_function_component[i][d] &&
2750 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2751 for (unsigned int qb = 0, q = 0; q < n_points;
2752 ++qb, q += n_lanes_user_interface)
2753 for (unsigned int v = 0;
2754 v < n_lanes_user_interface && q + v < n_points;
2755 ++v)
2756 ETT::access(this->gradients[qb],
2757 v,
2758 d,
2759 this->fe_values->shape_grad(i, q + v) * value);
2760 else if (this->nonzero_shape_function_component[i][d])
2761 for (unsigned int qb = 0, q = 0; q < n_points;
2762 ++qb, q += n_lanes_user_interface)
2763 for (unsigned int v = 0;
2764 v < n_lanes_user_interface && q + v < n_points;
2765 ++v)
2766 ETT::access(
2767 this->gradients[qb],
2768 v,
2769 d,
2770 this->fe_values->shape_grad_component(i, q + v, d) * value);
2771 }
2772 }
2773}
2774
2775
2776
2777template <int n_components_, int dim, int spacedim, typename Number>
2778template <bool is_linear>
2779inline void
2781 const EvaluationFlags::EvaluationFlags &integration_flags,
2782 const unsigned int n_shapes,
2783 const unsigned int qb,
2784 const vectorized_value_type value,
2786 vectorized_value_type *solution_values_vectorized_linear)
2787{
2788 if (integration_flags & EvaluationFlags::gradients)
2790 is_linear,
2791 dim,
2793 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2794 n_shapes,
2795 &value,
2796 gradient,
2797 is_linear ?
2798 solution_values_vectorized_linear :
2799 this->solution_renumbered_vectorized.data(),
2800 this->unit_point_ptr[qb],
2801 qb != 0);
2802 else
2804 dim,
2807 this->shapes.data() + qb * n_shapes,
2808 n_shapes,
2809 value,
2810 is_linear ? solution_values_vectorized_linear :
2811 this->solution_renumbered_vectorized.data(),
2812 this->unit_point_ptr[qb],
2813 qb != 0);
2814}
2815
2816
2817
2818template <int n_components_, int dim, int spacedim, typename Number>
2819template <bool is_linear, std::size_t stride_view>
2820inline void
2822 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2823 vectorized_value_type *solution_values_vectorized_linear,
2824 const bool sum_into_values)
2825{
2826 if (!sum_into_values && this->fe->n_components() > n_components)
2827 for (unsigned int i = 0; i < solution_values.size(); ++i)
2828 solution_values[i] = 0;
2829
2830 const unsigned int dofs_per_comp =
2831 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2832
2833 for (unsigned int comp = 0; comp < n_components; ++comp)
2834 {
2835 const std::size_t offset =
2836 (this->component_in_base_element + comp) * dofs_per_comp;
2837
2838 if (is_linear || this->renumber.empty())
2839 {
2840 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2841 if (sum_into_values)
2842 solution_values[i + offset] +=
2843 ETT::sum_value(comp,
2844 is_linear ?
2845 *(solution_values_vectorized_linear + i) :
2846 this->solution_renumbered_vectorized[i]);
2847 else
2848 solution_values[i + offset] =
2849 ETT::sum_value(comp,
2850 is_linear ?
2851 *(solution_values_vectorized_linear + i) :
2852 this->solution_renumbered_vectorized[i]);
2853 }
2854 else
2855 {
2856 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2857 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2858 if (sum_into_values)
2859 solution_values[renumber_ptr[i]] +=
2860 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2861 else
2862 solution_values[renumber_ptr[i]] =
2863 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2864 }
2865 }
2866}
2867
2868
2869
2870template <int n_components_, int dim, int spacedim, typename Number>
2871template <bool do_JxW, bool is_linear, std::size_t stride_view>
2872inline void
2874 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2875 const EvaluationFlags::EvaluationFlags &integration_flags,
2876 const bool sum_into_values)
2877{
2878 // zero out lanes of incomplete last quadrature point batch
2879 if constexpr (stride == 1)
2880 if (const unsigned int n_filled_lanes =
2881 this->n_q_points_scalar & (n_lanes_internal - 1);
2882 n_filled_lanes > 0)
2883 {
2884 if (integration_flags & EvaluationFlags::values)
2885 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2886 ETT::set_zero_value(this->values.back(), v);
2887 if (integration_flags & EvaluationFlags::gradients)
2888 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2889 ETT::set_zero_gradient(this->gradients.back(), v);
2890 }
2891
2892 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
2893 solution_values_vectorized_linear = {};
2894
2895 // loop over quadrature batches qb
2896 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2897
2898 const bool cartesian_cell =
2900 const bool affine_cell =
2902 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2903 {
2904 vectorized_value_type value = {};
2906
2907 if (integration_flags & EvaluationFlags::values)
2908 for (unsigned int v = 0, offset = qb * stride;
2909 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2910 ++v, ++offset)
2911 ETT::get_value(value,
2912 v,
2913 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
2914 this->values[offset]);
2915
2916 if (integration_flags & EvaluationFlags::gradients)
2917 for (unsigned int v = 0, offset = qb * stride;
2918 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2919 ++v, ++offset)
2920 {
2921 const gradient_type grad_w =
2922 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
2923 this->gradients[offset];
2924 ETT::get_gradient(
2925 gradient,
2926 v,
2927 cartesian_cell ?
2928 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
2929 grad_w) :
2931 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
2932 grad_w));
2933 }
2934
2935 compute_integrate_fast<is_linear>(
2936 integration_flags,
2937 n_shapes,
2938 qb,
2939 value,
2940 gradient,
2941 solution_values_vectorized_linear.data());
2942 }
2943
2944 // add between the lanes and write into the result
2945 finish_integrate_fast<is_linear>(solution_values,
2946 solution_values_vectorized_linear.data(),
2947 sum_into_values);
2948}
2949
2950
2951
2952template <int n_components_, int dim, int spacedim, typename Number>
2953template <bool do_JxW, std::size_t stride_view>
2954inline void
2956 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2957 const EvaluationFlags::EvaluationFlags &integration_flags,
2958 const bool sum_into_values)
2959{
2960 // slow path with FEValues
2961 Assert(this->fe_values.get() != nullptr,
2962 ExcMessage(
2963 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2964 if (!sum_into_values)
2965 for (unsigned int i = 0; i < solution_values.size(); ++i)
2966 solution_values[i] = 0;
2967
2968 const std::size_t n_points = this->fe_values->get_quadrature().size();
2969
2970 if (integration_flags & EvaluationFlags::values)
2971 {
2972 AssertIndexRange(this->n_q_points, this->values.size() + 1);
2973 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2974 {
2975 for (unsigned int d = 0; d < n_components; ++d)
2976 if (this->nonzero_shape_function_component[i][d] &&
2977 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2978 for (unsigned int qb = 0, q = 0; q < n_points;
2979 ++qb, q += n_lanes_user_interface)
2980 for (unsigned int v = 0;
2981 v < n_lanes_user_interface && q + v < n_points;
2982 ++v)
2983 solution_values[i] +=
2984 this->fe_values->shape_value(i, q + v) *
2985 ETT::access(this->values[qb], v, d) *
2986 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2987 else if (this->nonzero_shape_function_component[i][d])
2988 for (unsigned int qb = 0, q = 0; q < n_points;
2989 ++qb, q += n_lanes_user_interface)
2990 for (unsigned int v = 0;
2991 v < n_lanes_user_interface && q + v < n_points;
2992 ++v)
2993 solution_values[i] +=
2994 this->fe_values->shape_value_component(i, q + v, d) *
2995 ETT::access(this->values[qb], v, d) *
2996 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2997 }
2998 }
2999
3000 if (integration_flags & EvaluationFlags::gradients)
3001 {
3002 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
3003 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3004 {
3005 for (unsigned int d = 0; d < n_components; ++d)
3006 if (this->nonzero_shape_function_component[i][d] &&
3007 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3008 for (unsigned int qb = 0, q = 0; q < n_points;
3009 ++qb, q += n_lanes_user_interface)
3010 for (unsigned int v = 0;
3011 v < n_lanes_user_interface && q + v < n_points;
3012 ++v)
3013 solution_values[i] +=
3014 this->fe_values->shape_grad(i, q + v) *
3015 ETT::access(this->gradients[qb], v, d) *
3016 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3017 else if (this->nonzero_shape_function_component[i][d])
3018 for (unsigned int qb = 0, q = 0; q < n_points;
3019 ++qb, q += n_lanes_user_interface)
3020 for (unsigned int v = 0;
3021 v < n_lanes_user_interface && q + v < n_points;
3022 ++v)
3023 solution_values[i] +=
3024 this->fe_values->shape_grad_component(i, q + v, d) *
3025 ETT::access(this->gradients[qb], v, d) *
3026 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3027 }
3028 }
3029}
3030
3031
3032
3033template <int n_components_, int dim, int spacedim, typename Number>
3034template <bool do_JxW, std::size_t stride_view>
3035void
3037 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3038 const EvaluationFlags::EvaluationFlags &integration_flags,
3039 const bool sum_into_values)
3040{
3041 if (!this->is_reinitialized)
3042 reinit();
3043
3044 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3045
3046 if (this->n_q_points == 0 || // no evaluation points provided
3047 !((integration_flags & EvaluationFlags::values) ||
3048 (integration_flags &
3049 EvaluationFlags::gradients))) // no integration flags
3050 {
3051 if (!sum_into_values)
3052 for (unsigned int i = 0; i < solution_values.size(); ++i)
3053 solution_values[i] = 0;
3054 return;
3055 }
3056
3057 Assert(
3058 !do_JxW || this->JxW_ptr != nullptr,
3059 ExcMessage(
3060 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3061
3062 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3063 if (this->fast_path)
3064 {
3065 if (this->use_linear_path)
3066 integrate_fast<do_JxW, true>(solution_values,
3067 integration_flags,
3068 sum_into_values);
3069 else
3070 integrate_fast<do_JxW, false>(solution_values,
3071 integration_flags,
3072 sum_into_values);
3073 }
3074 else
3075 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3076}
3077
3078
3079
3080template <int n_components_, int dim, int spacedim, typename Number>
3083 const unsigned int point_index) const
3084{
3085 AssertIndexRange(point_index, this->n_q_points);
3086 Assert(this->normal_ptr != nullptr,
3088 ExcFEPointEvaluationAccessToUninitializedMappingField(
3089 "update_normal_vectors"));
3090 if (this->is_interior)
3091 return this->normal_ptr[point_index];
3092 else
3093 return -this->normal_ptr[point_index];
3094}
3095
3096
3097
3098template <int n_components_, int dim, int spacedim, typename Number>
3103 const bool is_interior,
3104 const unsigned int first_selected_component)
3105 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3106 mapping_info,
3107 fe,
3108 first_selected_component,
3109 is_interior)
3110{}
3111
3112
3113
3114template <int n_components_, int dim, int spacedim, typename Number>
3115inline void
3117 const unsigned int cell_index,
3118 const unsigned int face_number)
3119{
3120 this->current_cell_index = cell_index;
3121 this->current_face_number = face_number;
3122
3123 if (this->use_linear_path)
3124 this->template do_reinit<true, true>();
3125 else
3126 this->template do_reinit<true, false>();
3127}
3128
3129
3130
3131template <int n_components_, int dim, int spacedim, typename Number>
3132inline void
3134 const unsigned int face_index)
3135{
3136 this->current_cell_index = face_index;
3137 this->current_face_number =
3138 this->mapping_info->get_face_number(face_index, this->is_interior);
3139
3140 if (this->use_linear_path)
3141 this->template do_reinit<true, true>();
3142 else
3143 this->template do_reinit<true, false>();
3144}
3145
3146
3147
3148template <int n_components_, int dim, int spacedim, typename Number>
3149template <std::size_t stride_view>
3150void
3153 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3154{
3155 Assert(this->is_reinitialized, ExcMessage("Is not reinitialized!"));
3156
3157 if (this->n_q_points == 0)
3158 return;
3159
3160 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3161
3162 if (!((evaluation_flags & EvaluationFlags::values) ||
3163 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3164 return;
3165
3166 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3167
3168 if (this->use_linear_path)
3169 do_evaluate<true>(solution_values, evaluation_flags);
3170 else
3171 do_evaluate<false>(solution_values, evaluation_flags);
3172}
3173
3174
3175
3176template <int n_components_, int dim, int spacedim, typename Number>
3177void
3179 const ArrayView<const ScalarNumber> &solution_values,
3180 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3181{
3182 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3183 solution_values.size()),
3184 evaluation_flags);
3185}
3186
3187
3188
3189template <int n_components_, int dim, int spacedim, typename Number>
3190template <bool is_linear, std::size_t stride_view>
3191void
3194 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3195{
3196 const unsigned int dofs_per_comp =
3197 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3198
3199 const ScalarNumber *input;
3200 if (stride_view == 1 && this->component_in_base_element == 0 &&
3201 (is_linear || this->renumber.empty()))
3202 input = solution_values.data();
3203 else
3204 {
3205 for (unsigned int comp = 0; comp < n_components; ++comp)
3206 {
3207 const std::size_t offset =
3208 (this->component_in_base_element + comp) * dofs_per_comp;
3209
3210 if (is_linear || this->renumber.empty())
3211 {
3212 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3213 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3214 solution_values[i + offset];
3215 }
3216 else
3217 {
3218 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3219 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3220 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3221 solution_values[renumber_ptr[i]];
3222 }
3223 }
3224 input = this->scratch_data_scalar.data();
3225 }
3226
3227 ScalarNumber *output =
3228 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3229
3230 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3231 template interpolate<true, false>(n_components,
3232 evaluation_flags,
3233 this->shape_info,
3234 input,
3235 output,
3236 this->current_face_number);
3237
3238 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3239}
3240
3241
3242
3243template <int n_components_, int dim, int spacedim, typename Number>
3244template <std::size_t stride_view>
3245void
3247 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3248 const EvaluationFlags::EvaluationFlags &integration_flags,
3249 const bool sum_into_values)
3250{
3251 Assert(this->is_reinitialized, ExcMessage("Is not reinitialized!"));
3252
3253 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3254
3255 if (this->n_q_points == 0 || // no evaluation points provided
3256 !((integration_flags & EvaluationFlags::values) ||
3257 (integration_flags &
3258 EvaluationFlags::gradients))) // no integration flags
3259 {
3260 if (!sum_into_values)
3261 for (unsigned int i = 0; i < solution_values.size(); ++i)
3262 solution_values[i] = 0;
3263 return;
3264 }
3265
3266 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3267
3268 if (this->use_linear_path)
3269 do_integrate<true, true>(solution_values,
3270 integration_flags,
3271 sum_into_values);
3272 else
3273 do_integrate<true, false>(solution_values,
3274 integration_flags,
3275 sum_into_values);
3276}
3277
3278
3279
3280template <int n_components_, int dim, int spacedim, typename Number>
3281void
3283 const ArrayView<ScalarNumber> &solution_values,
3284 const EvaluationFlags::EvaluationFlags &integration_flags,
3285 const bool sum_into_values)
3286{
3287 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3288 solution_values.size()),
3289 integration_flags,
3290 sum_into_values);
3291}
3292
3293
3294
3295template <int n_components_, int dim, int spacedim, typename Number>
3296template <std::size_t stride_view>
3297void
3299 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3300 const EvaluationFlags::EvaluationFlags &integration_flags,
3301 const bool sum_into_values)
3302{
3303 Assert(this->is_reinitialized, ExcMessage("Is not reinitialized!"));
3304
3305 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3306
3307 if (this->n_q_points == 0 || // no evaluation points provided
3308 !((integration_flags & EvaluationFlags::values) ||
3309 (integration_flags &
3310 EvaluationFlags::gradients))) // no integration flags
3311 {
3312 if (!sum_into_values)
3313 for (unsigned int i = 0; i < solution_values.size(); ++i)
3314 solution_values[i] = 0;
3315 return;
3316 }
3317
3318 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3319
3320 if (this->use_linear_path)
3321 do_integrate<false, true>(solution_values,
3322 integration_flags,
3323 sum_into_values);
3324 else
3325 do_integrate<false, false>(solution_values,
3326 integration_flags,
3327 sum_into_values);
3328}
3329
3330
3331
3332template <int n_components_, int dim, int spacedim, typename Number>
3333void
3335 const ArrayView<ScalarNumber> &solution_values,
3336 const EvaluationFlags::EvaluationFlags &integration_flags,
3337 const bool sum_into_values)
3338{
3339 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3340 solution_values.size()),
3341 integration_flags,
3342 sum_into_values);
3343}
3344
3345
3346
3347template <int n_components_, int dim, int spacedim, typename Number>
3348template <bool do_JxW, bool is_linear, std::size_t stride_view>
3349void
3351 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3352 const EvaluationFlags::EvaluationFlags &integration_flags,
3353 const bool sum_into_values)
3354{
3355 if (!sum_into_values && this->fe->n_components() > n_components)
3356 for (unsigned int i = 0; i < solution_values.size(); ++i)
3357 solution_values[i] = 0;
3358
3359 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3360 integration_flags,
3361 false);
3362
3363 ScalarNumber *input = this->scratch_data_scalar.begin();
3364
3365 if (stride_view == 1 && this->component_in_base_element == 0 &&
3366 (is_linear || this->renumber.empty()))
3367 {
3368 if (sum_into_values)
3369 internal::
3370 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3371 template interpolate<false, true>(n_components,
3372 integration_flags,
3373 this->shape_info,
3374 input,
3375 solution_values.data(),
3376 this->current_face_number);
3377 else
3378 internal::
3379 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3380 template interpolate<false, false>(n_components,
3381 integration_flags,
3382 this->shape_info,
3383 input,
3384 solution_values.data(),
3385 this->current_face_number);
3386 }
3387 else
3388 {
3389 const unsigned int dofs_per_comp_face =
3390 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3391
3392 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3393 ScalarNumber *output = input + size_input;
3394
3395 internal::
3396 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3397 template interpolate<false, false>(n_components,
3398 integration_flags,
3399 this->shape_info,
3400 input,
3401 output,
3402 this->current_face_number);
3403
3404 const unsigned int dofs_per_comp =
3405 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3406
3407 for (unsigned int comp = 0; comp < n_components; ++comp)
3408 {
3409 const std::size_t offset =
3410 (this->component_in_base_element + comp) * dofs_per_comp;
3411
3412 if (is_linear || this->renumber.empty())
3413 {
3414 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3415 if (sum_into_values)
3416 solution_values[i + offset] +=
3417 output[i + comp * dofs_per_comp];
3418 else
3419 solution_values[i + offset] =
3420 output[i + comp * dofs_per_comp];
3421 }
3422 else
3423 {
3424 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3425 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3426 if (sum_into_values)
3427 solution_values[renumber_ptr[i]] +=
3428 output[i + comp * dofs_per_comp];
3429 else
3430 solution_values[renumber_ptr[i]] =
3431 output[i + comp * dofs_per_comp];
3432 }
3433 }
3434 }
3435}
3436
3437
3438
3439template <int n_components_, int dim, int spacedim, typename Number>
3440template <int stride_face_dof>
3441void
3443 const ScalarNumber *face_dof_values,
3444 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3445{
3446 if (this->use_linear_path)
3447 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3448 evaluation_flags);
3449 else
3450 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3451 evaluation_flags);
3452}
3453
3454
3455
3456template <int n_components_, int dim, int spacedim, typename Number>
3457template <bool is_linear, int stride_face_dof>
3458inline void
3460 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3461 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3462{
3463 const scalar_value_type *face_dof_values_ptr;
3464 if constexpr (n_components == 1)
3465 face_dof_values_ptr = face_dof_values;
3466 else
3467 {
3468 const unsigned int dofs_per_comp_face =
3469 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3470 for (unsigned int comp = 0; comp < n_components; ++comp)
3471 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3472 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3473 stride_face_dof],
3474 comp,
3475 this->solution_renumbered[i]);
3476
3477 face_dof_values_ptr = this->solution_renumbered.data();
3478 }
3479
3480 constexpr int stride_face_dof_actual =
3481 n_components == 1 ? stride_face_dof : 1;
3482
3483 // loop over quadrature batches qb
3484 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3485
3486 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3487 {
3490
3491 if (evaluation_flags & EvaluationFlags::gradients)
3492 {
3493 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3494 is_linear ?
3496 dim - 1,
3499 2,
3500 stride_face_dof_actual>(face_dof_values_ptr,
3501 this->unit_point_faces_ptr[qb]) :
3503 dim - 1,
3506 2,
3507 false,
3508 stride_face_dof_actual>(this->shapes_faces.data() +
3509 qb * n_shapes,
3510 n_shapes,
3511 face_dof_values_ptr);
3512
3513 value = interpolated_value[dim - 1];
3514 // reorder derivative from tangential/normal derivatives into tensor
3515 // in physical coordinates
3516 if (this->current_face_number / 2 == 0)
3517 {
3518 gradient[0] = interpolated_value[dim];
3519 if (dim > 1)
3520 gradient[1] = interpolated_value[0];
3521 if (dim > 2)
3522 gradient[2] = interpolated_value[1];
3523 }
3524 else if (this->current_face_number / 2 == 1)
3525 {
3526 if (dim > 1)
3527 gradient[1] = interpolated_value[dim];
3528 if (dim == 3)
3529 {
3530 gradient[0] = interpolated_value[1];
3531 gradient[2] = interpolated_value[0];
3532 }
3533 else if (dim == 2)
3534 gradient[0] = interpolated_value[0];
3535 else
3537 }
3538 else if (this->current_face_number / 2 == 2)
3539 {
3540 if (dim > 2)
3541 {
3542 gradient[0] = interpolated_value[0];
3543 gradient[1] = interpolated_value[1];
3544 gradient[2] = interpolated_value[dim];
3545 }
3546 else
3548 }
3549 else
3551 }
3552 else
3553 {
3554 value = is_linear ?
3556 dim - 1,
3559 stride_face_dof_actual>(face_dof_values_ptr,
3560 this->unit_point_faces_ptr[qb]) :
3562 dim - 1,
3565 false,
3566 stride_face_dof_actual>(this->shapes_faces.data() +
3567 qb * n_shapes,
3568 n_shapes,
3569 face_dof_values_ptr);
3570 }
3571
3572 if (evaluation_flags & EvaluationFlags::values)
3573 {
3574 for (unsigned int v = 0, offset = qb * stride;
3575 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3576 ++v, ++offset)
3577 ETT::set_value(value, v, this->values[offset]);
3578 }
3579 if (evaluation_flags & EvaluationFlags::gradients)
3580 {
3581 Assert(this->update_flags & update_gradients ||
3582 this->update_flags & update_inverse_jacobians,
3584
3585 for (unsigned int v = 0, offset = qb * stride;
3586 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3587 ++v, ++offset)
3588 {
3589 unit_gradient_type unit_gradient;
3590 ETT::set_gradient(gradient, v, unit_gradient);
3591 this->gradients[offset] =
3592 this->cell_type <=
3595 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3597 this
3598 ->inverse_jacobian_ptr[this->cell_type <=
3600 GeometryType::affine ?
3601 0 :
3602 offset]
3603 .transpose(),
3604 unit_gradient);
3605 }
3606 }
3607 }
3608}
3609
3610
3611
3612template <int n_components_, int dim, int spacedim, typename Number>
3613template <int stride_face_dof>
3614void
3616 ScalarNumber *face_dof_values,
3617 const EvaluationFlags::EvaluationFlags &integration_flags,
3618 const bool sum_into_values)
3619{
3620 if (this->use_linear_path)
3621 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3622 integration_flags,
3623 sum_into_values);
3624 else
3625 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3626 integration_flags,
3627 sum_into_values);
3628}
3629
3630
3631
3632template <int n_components_, int dim, int spacedim, typename Number>
3633template <bool do_JxW, bool is_linear, int stride_face_dof>
3634inline void
3637 ScalarNumber *face_dof_values,
3638 const EvaluationFlags::EvaluationFlags &integration_flags,
3639 const bool sum_into_values)
3640{
3641 // zero out lanes of incomplete last quadrature point batch
3642 if constexpr (stride == 1)
3643 if (const unsigned int n_filled_lanes =
3644 this->n_q_points_scalar & (n_lanes_internal - 1);
3645 n_filled_lanes > 0)
3646 {
3647 if (integration_flags & EvaluationFlags::values)
3648 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3649 ETT::set_zero_value(this->values.back(), v);
3650 if (integration_flags & EvaluationFlags::gradients)
3651 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3652 ETT::set_zero_gradient(this->gradients.back(), v);
3653 }
3654
3655 std::array<vectorized_value_type,
3656 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3657 solution_values_vectorized_linear = {};
3658
3659 // loop over quadrature batches qb
3660 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3661
3662 const bool cartesian_cell =
3664 const bool affine_cell =
3666 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3667 {
3668 vectorized_value_type value = {};
3670
3671 if (integration_flags & EvaluationFlags::values)
3672 for (unsigned int v = 0, offset = qb * stride;
3673 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3674 ++v, ++offset)
3675 ETT::get_value(value,
3676 v,
3677 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3678 this->values[offset]);
3679
3680 if (integration_flags & EvaluationFlags::gradients)
3681 for (unsigned int v = 0, offset = qb * stride;
3682 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3683 ++v, ++offset)
3684 {
3685 const auto grad_w =
3686 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3687 this->gradients[offset];
3688 ETT::get_gradient(
3689 gradient,
3690 v,
3691 cartesian_cell ?
3692 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3693 grad_w) :
3695 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3696 grad_w));
3697 }
3698
3699 if (integration_flags & EvaluationFlags::gradients)
3700 {
3701 std::array<vectorized_value_type, 2> value_face = {};
3702 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3703
3704 value_face[0] = value;
3705 // fill derivative in physical coordinates into tangential/normal
3706 // derivatives
3707 if (this->current_face_number / 2 == 0)
3708 {
3709 value_face[1] = gradient[0];
3710 if (dim > 1)
3711 gradient_in_face[0] = gradient[1];
3712 if (dim > 2)
3713 gradient_in_face[1] = gradient[2];
3714 }
3715 else if (this->current_face_number / 2 == 1)
3716 {
3717 if (dim > 1)
3718 value_face[1] = gradient[1];
3719 if (dim == 3)
3720 {
3721 gradient_in_face[0] = gradient[2];
3722 gradient_in_face[1] = gradient[0];
3723 }
3724 else if (dim == 2)
3725 gradient_in_face[0] = gradient[0];
3726 else
3728 }
3729 else if (this->current_face_number / 2 == 2)
3730 {
3731 if (dim > 2)
3732 {
3733 value_face[1] = gradient[2];
3734 gradient_in_face[0] = gradient[0];
3735 gradient_in_face[1] = gradient[1];
3736 }
3737 else
3739 }
3740 else
3742
3744 is_linear,
3745 dim - 1,
3748 2>(this->shapes_faces.data() + qb * n_shapes,
3749 n_shapes,
3750 value_face.data(),
3751 gradient_in_face,
3752 is_linear ? solution_values_vectorized_linear.data() :
3753 this->solution_renumbered_vectorized.data(),
3754 this->unit_point_faces_ptr[qb],
3755 qb != 0);
3756 }
3757 else
3759 dim - 1,
3762 this->shapes_faces.data() + qb * n_shapes,
3763 n_shapes,
3764 value,
3765 is_linear ? solution_values_vectorized_linear.data() :
3766 this->solution_renumbered_vectorized.data(),
3767 this->unit_point_faces_ptr[qb],
3768 qb != 0);
3769 }
3770
3771 const unsigned int dofs_per_comp_face =
3772 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3773
3774 for (unsigned int comp = 0; comp < n_components; ++comp)
3775 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3776 if (sum_into_values)
3777 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3778 stride_face_dof] +=
3779 ETT::sum_value(comp,
3780 is_linear ?
3781 *(solution_values_vectorized_linear.data() + i) :
3782 this->solution_renumbered_vectorized[i]);
3783 else
3784 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3785 ETT::sum_value(comp,
3786 is_linear ?
3787 *(solution_values_vectorized_linear.data() + i) :
3788 this->solution_renumbered_vectorized[i]);
3789}
3790
3791
3792
3793template <int n_components_, int dim, int spacedim, typename Number>
3796 const unsigned int point_index) const
3797{
3798 AssertIndexRange(point_index, this->n_q_points);
3799 Assert(this->normal_ptr != nullptr,
3801 ExcFEPointEvaluationAccessToUninitializedMappingField(
3802 "update_normal_vectors"));
3803 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3804 {
3806 for (unsigned int d = 0; d < dim; ++d)
3807 normal[d] =
3808 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3809 0);
3810 if (this->is_interior)
3811 return normal;
3812 else
3813 return -normal;
3814 }
3815 else
3816 {
3817 if (this->is_interior)
3818 return this->normal_ptr[point_index];
3819 else
3820 return -(this->normal_ptr[point_index]);
3821 }
3822}
3823
3825
3826#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
boost::signals2::connection connection_is_reinitialized
const Tensor< 1, spacedim, Number > * normal_ptr
DerivativeForm< 1, spacedim, dim, Number > inverse_jacobian(const unsigned int point_index) const
AlignedVector< vectorized_value_type > solution_renumbered_vectorized
static constexpr std::size_t n_lanes_internal
const DerivativeForm< 1, spacedim, dim, Number > * inverse_jacobian_ptr
AlignedVector<::ndarray< VectorizedArrayType, 2, dim > > shapes
typename ETT::value_type value_type
scalar_value_type integrate_value() const
void submit_gradient(const gradient_type &, const unsigned int point_index)
void setup(const unsigned int first_selected_component)
FEPointEvaluationBase(NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const unsigned int first_selected_component=0, const bool is_interior=true)
typename ETT::scalar_value_type scalar_value_type
static constexpr unsigned int dimension
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::vectorized_value_type vectorized_value_type
static constexpr unsigned int n_components
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
FEPointEvaluationBase(FEPointEvaluationBase &other) noexcept
DerivativeForm< 1, dim, spacedim, Number > jacobian(const unsigned int point_index) const
typename ETT::real_gradient_type gradient_type
const Point< spacedim, Number > * real_point_ptr
unsigned int component_in_base_element
SmartPointer< const Mapping< dim, spacedim > > mapping
internal::MatrixFreeFunctions::ShapeInfo< ScalarNumber > shape_info
void submit_value(const value_type &value, const unsigned int point_index)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static constexpr std::size_t stride
void integrate_slow(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr unsigned int dimension
typename ETT::value_type value_type
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename ETT::scalar_value_type scalar_value_type
typename ETT::unit_gradient_type unit_gradient_type
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void prepare_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values)
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
typename ETT::vectorized_value_type vectorized_value_type
void compute_integrate_fast(const EvaluationFlags::EvaluationFlags &integration_flags, const unsigned int n_shapes, const unsigned int qb, const vectorized_value_type value, const interface_vectorized_unit_gradient_type gradient, vectorized_value_type *solution_values_vectorized_linear)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::real_gradient_type gradient_type
FEPointEvaluation(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void evaluate_slow(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr std::size_t n_lanes_internal
static constexpr unsigned int n_components
void finish_integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, vectorized_value_type *solution_values_vectorized_linear, const bool sum_into_values)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_unit_gradient_type &gradient)
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr std::size_t n_lanes_user_interface
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
value_type * data() const noexcept
Definition array_view.h:846
std::size_t size() const
Definition array_view.h:867
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_diagonal_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcFEPointEvaluationAccessToUninitializedMappingField(std::string arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_ASSERT_UNREACHABLE()
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool do_add)
void compute_values_of_array_in_pairs(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p0, const Point< dim, Number > &p1)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const Number *values, const Point< dim, Number2 > &p)
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool do_add)
static const unsigned int invalid_unsigned_int
Definition types.h:220
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
STL namespace.
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
static void set_zero_gradient(unit_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static Tensor< 1, dim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, dim, ScalarNumber > &shape_gradient)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int, const vectorized_unit_gradient_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int vector_lane, scalar_unit_gradient_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static scalar_value_type sum_value(const vectorized_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int, vectorized_unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int, const ScalarNumber &shape_value)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int vector_lane, const scalar_unit_gradient_type &result)
static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const DerivativeForm< 1, dim, n_components, Number > &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static scalar_value_type sum_value(const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
Tensor< 1, n_components, ScalarNumber > scalar_value_type
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
Tensor< 1, n_components, Tensor< 1, dim, VectorizedArrayType > > vectorized_unit_gradient_type
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
std::conditional_t< n_components==spacedim, Tensor< 2, spacedim, Number >, Tensor< 1, n_components, Tensor< 1, spacedim, Number > > > real_gradient_type
Tensor< 1, n_components, VectorizedArrayType > vectorized_value_type
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
Tensor< 1, n_components, Tensor< 1, dim, Number > > unit_gradient_type
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static constexpr std::size_t width()
static constexpr std::size_t stride()
static value_type & get(value_type &value, unsigned int c)