Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00: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
819protected:
820 static constexpr std::size_t n_lanes_user_interface =
822 static constexpr std::size_t n_lanes_internal =
824 static constexpr std::size_t stride =
826
835 void
836 setup(const unsigned int first_selected_component);
837
843 template <bool is_face, bool is_linear>
844 void
846
850 const unsigned int n_q_batches;
851
855 const unsigned int n_q_points;
856
860 const unsigned int n_q_points_scalar;
861
866
871
876 std::vector<Polynomials::Polynomial<double>> poly;
877
882
887 std::vector<unsigned int> renumber;
888
896 std::vector<scalar_value_type> solution_renumbered;
897
905
910
914 std::vector<value_type> values;
915
919 std::vector<gradient_type> gradients;
920
926
932
938
944
950
956
962 const Number *JxW_ptr;
963
968
973 unsigned int dofs_per_component;
974
981
986
991
997 std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
998
1003
1007 std::shared_ptr<FEValues<dim, spacedim>> fe_values;
1008
1012 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1014
1020
1025
1030
1035
1040 boost::signals2::connection connection_is_reinitialized;
1041
1046
1052
1058
1059 const bool is_interior;
1060};
1061
1062
1063
1094template <int n_components_,
1095 int dim,
1096 int spacedim = dim,
1097 typename Number = double>
1099 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1100{
1101public:
1102 static constexpr unsigned int dimension = dim;
1103 static constexpr unsigned int n_components = n_components_;
1104
1105 using number_type = Number;
1106
1109 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1111 using ETT = typename internal::FEPointEvaluation::
1112 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1113 using value_type = typename ETT::value_type;
1114 using scalar_value_type = typename ETT::scalar_value_type;
1115 using vectorized_value_type = typename ETT::vectorized_value_type;
1116 using unit_gradient_type = typename ETT::unit_gradient_type;
1117 using gradient_type = typename ETT::real_gradient_type;
1119 typename ETT::interface_vectorized_unit_gradient_type;
1120
1142 const unsigned int first_selected_component = 0);
1143
1163 const unsigned int first_selected_component = 0);
1164
1176 void
1178 const ArrayView<const Point<dim>> &unit_points);
1179
1184 void
1185 reinit();
1186
1191 void
1192 reinit(const unsigned int cell_index);
1193
1194
1206 template <std::size_t stride_view>
1207 void
1208 evaluate(
1210 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1211
1223 void
1224 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1225 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1226
1249 template <std::size_t stride_view>
1250 void
1252 const EvaluationFlags::EvaluationFlags &integration_flags,
1253 const bool sum_into_values = false);
1254
1277 void
1278 integrate(const ArrayView<ScalarNumber> &solution_values,
1279 const EvaluationFlags::EvaluationFlags &integration_flags,
1280 const bool sum_into_values = false);
1281
1308 template <std::size_t stride_view>
1309 void
1311 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1312 const EvaluationFlags::EvaluationFlags &integration_flags,
1313 const bool sum_into_values = false);
1314
1341 void
1342 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1343 const EvaluationFlags::EvaluationFlags &integration_flags,
1344 const bool sum_into_values = false);
1345
1352 normal_vector(const unsigned int point_index) const;
1353
1354private:
1355 static constexpr std::size_t n_lanes_user_interface =
1357 static constexpr std::size_t n_lanes_internal =
1359 static constexpr std::size_t stride =
1361
1366 template <bool is_linear, std::size_t stride_view>
1367 void
1370
1375 template <bool is_linear, std::size_t stride_view>
1376 void
1379 const EvaluationFlags::EvaluationFlags &evaluation_flags,
1380 const unsigned int n_shapes,
1381 const unsigned int qb,
1382 vectorized_value_type &value,
1384
1388 template <bool is_linear, std::size_t stride_view>
1389 void
1392 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1393
1397 template <std::size_t stride_view>
1398 void
1401 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1402
1408 template <bool is_linear>
1409 void
1411 const EvaluationFlags::EvaluationFlags &integration_flags,
1412 const unsigned int n_shapes,
1413 const unsigned int qb,
1414 const vectorized_value_type value,
1416 vectorized_value_type *solution_values_vectorized_linear);
1417
1423 template <bool is_linear, std::size_t stride_view>
1424 void
1426 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1427 vectorized_value_type *solution_values_vectorized_linear,
1428 const bool sum_into_values);
1429
1433 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1434 void
1436 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1437 const EvaluationFlags::EvaluationFlags &integration_flags,
1438 const bool sum_into_values);
1439
1443 template <bool do_JxW, std::size_t stride_view>
1444 void
1446 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1447 const EvaluationFlags::EvaluationFlags &integration_flags,
1448 const bool sum_into_values);
1449
1453 template <bool do_JxW, std::size_t stride_view>
1454 void
1456 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1457 const EvaluationFlags::EvaluationFlags &integration_flags,
1458 const bool sum_into_values);
1459};
1460
1461
1462
1480template <int n_components_,
1481 int dim,
1482 int spacedim = dim,
1483 typename Number = double>
1485 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1486{
1487public:
1488 static constexpr unsigned int dimension = dim;
1489 static constexpr unsigned int n_components = n_components_;
1490
1491 using number_type = Number;
1492
1495 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1497 using ETT = typename internal::FEPointEvaluation::
1498 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1499 using value_type = typename ETT::value_type;
1500 using scalar_value_type = typename ETT::scalar_value_type;
1501 using vectorized_value_type = typename ETT::vectorized_value_type;
1502 using unit_gradient_type = typename ETT::unit_gradient_type;
1503 using gradient_type = typename ETT::real_gradient_type;
1505 typename ETT::interface_vectorized_unit_gradient_type;
1506
1513 const bool is_interior = true,
1514 const unsigned int first_selected_component = 0);
1515
1520 void
1521 reinit(const unsigned int cell_index, const unsigned int face_number);
1522
1527 void
1528 reinit(const unsigned int face_index);
1529
1541 template <std::size_t stride_view>
1542 void
1543 evaluate(
1545 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1546
1558 void
1559 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1560 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1561
1584 template <std::size_t stride_view>
1585 void
1587 const EvaluationFlags::EvaluationFlags &integration_flags,
1588 const bool sum_into_values = false);
1589
1612 void
1613 integrate(const ArrayView<ScalarNumber> &solution_values,
1614 const EvaluationFlags::EvaluationFlags &integration_flags,
1615 const bool sum_into_values = false);
1616
1639 template <std::size_t stride_view>
1640 void
1642 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1643 const EvaluationFlags::EvaluationFlags &integration_flags,
1644 const bool sum_into_values = false);
1645
1668 void
1669 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1670 const EvaluationFlags::EvaluationFlags &integration_flags,
1671 const bool sum_into_values = false);
1672
1679 template <int stride_face_dof = VectorizedArrayType::size()>
1680 void
1681 evaluate_in_face(const ScalarNumber *face_dof_values,
1682 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1683
1690 template <int stride_face_dof = VectorizedArrayType::size()>
1691 void
1692 integrate_in_face(ScalarNumber *face_dof_values,
1693 const EvaluationFlags::EvaluationFlags &integration_flags,
1694 const bool sum_into_values = false);
1695
1702 normal_vector(const unsigned int point_index) const;
1703
1704private:
1705 static constexpr std::size_t n_lanes_user_interface =
1707 static constexpr std::size_t n_lanes_internal =
1709 static constexpr std::size_t stride =
1711
1712 template <bool is_linear, std::size_t stride_view>
1713 void
1716 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1717
1718 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1719 void
1721 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1722 const EvaluationFlags::EvaluationFlags &integration_flags,
1723 const bool sum_into_values);
1724
1729 template <bool is_linear, int stride_face_dof>
1730 void
1731 do_evaluate_in_face(const ScalarNumber *face_dof_values,
1732 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1733
1738 template <bool do_JxW, bool is_linear, int stride_face_dof>
1739 void
1741 ScalarNumber *face_dof_values,
1742 const EvaluationFlags::EvaluationFlags &integration_flags,
1743 const bool sum_into_values);
1744};
1745
1746
1747
1748// ----------------------- template and inline function ----------------------
1749
1750
1751template <int n_components_, int dim, int spacedim, typename Number>
1755 const UpdateFlags update_flags,
1756 const unsigned int first_selected_component)
1757 : n_q_batches(numbers::invalid_unsigned_int)
1758 , n_q_points(numbers::invalid_unsigned_int)
1759 , n_q_points_scalar(numbers::invalid_unsigned_int)
1760 , mapping(&mapping)
1761 , fe(&fe)
1762 , JxW_ptr(nullptr)
1763 , update_flags(update_flags)
1764 , mapping_info_on_the_fly(
1765 std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1766 mapping,
1767 update_flags))
1768 , mapping_info(mapping_info_on_the_fly.get())
1769 , current_cell_index(numbers::invalid_unsigned_int)
1770 , current_face_number(numbers::invalid_unsigned_int)
1771 , is_reinitialized(false)
1772 , is_interior(true)
1773{
1774 setup(first_selected_component);
1775}
1776
1777
1778
1779template <int n_components_, int dim, int spacedim, typename Number>
1784 const unsigned int first_selected_component,
1785 const bool is_interior)
1786 : n_q_batches(numbers::invalid_unsigned_int)
1787 , n_q_points(numbers::invalid_unsigned_int)
1788 , n_q_points_scalar(numbers::invalid_unsigned_int)
1789 , mapping(&mapping_info.get_mapping())
1790 , fe(&fe)
1791 , JxW_ptr(nullptr)
1792 , update_flags(mapping_info.get_update_flags())
1793 , mapping_info(&mapping_info)
1794 , current_cell_index(numbers::invalid_unsigned_int)
1795 , current_face_number(numbers::invalid_unsigned_int)
1796 , is_reinitialized(false)
1797 , is_interior(is_interior)
1798{
1799 setup(first_selected_component);
1800 connection_is_reinitialized = mapping_info.connect_is_reinitialized(
1801 [this]() { this->is_reinitialized = false; });
1802}
1803
1804
1805
1806template <int n_components_, int dim, int spacedim, typename Number>
1810 : n_q_batches(other.n_q_batches)
1811 , n_q_points(other.n_q_points)
1812 , n_q_points_scalar(other.n_q_points_scalar)
1813 , mapping(other.mapping)
1814 , fe(other.fe)
1815 , poly(other.poly)
1816 , use_linear_path(other.use_linear_path)
1817 , renumber(other.renumber)
1818 , solution_renumbered(other.solution_renumbered)
1819 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1820 , values(other.values)
1821 , gradients(other.gradients)
1822 , dofs_per_component(other.dofs_per_component)
1823 , dofs_per_component_face(other.dofs_per_component_face)
1824 , component_in_base_element(other.component_in_base_element)
1825 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1826 , update_flags(other.update_flags)
1827 , fe_values(other.fe_values)
1828 , mapping_info_on_the_fly(
1829 other.mapping_info_on_the_fly ?
1831 *mapping,
1832 update_flags) :
1833 nullptr)
1834 , mapping_info(other.mapping_info)
1835 , current_cell_index(other.current_cell_index)
1836 , current_face_number(other.current_face_number)
1837 , fast_path(other.fast_path)
1838 , is_reinitialized(false)
1839 , shapes(other.shapes)
1840 , shapes_faces(other.shapes_faces)
1841 , is_interior(other.is_interior)
1842{
1843 connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1844 [this]() { this->is_reinitialized = false; });
1845}
1846
1847
1848
1849template <int n_components_, int dim, int spacedim, typename Number>
1852 : n_q_batches(other.n_q_batches)
1853 , n_q_points(other.n_q_points)
1854 , n_q_points_scalar(other.n_q_points_scalar)
1855 , mapping(other.mapping)
1856 , fe(other.fe)
1857 , poly(other.poly)
1858 , use_linear_path(other.use_linear_path)
1859 , renumber(other.renumber)
1860 , solution_renumbered(other.solution_renumbered)
1861 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1862 , values(other.values)
1863 , gradients(other.gradients)
1864 , dofs_per_component(other.dofs_per_component)
1865 , dofs_per_component_face(other.dofs_per_component_face)
1866 , component_in_base_element(other.component_in_base_element)
1867 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1868 , update_flags(other.update_flags)
1869 , fe_values(other.fe_values)
1870 , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1871 , mapping_info(other.mapping_info)
1872 , current_cell_index(other.current_cell_index)
1873 , current_face_number(other.current_face_number)
1874 , fast_path(other.fast_path)
1875 , is_reinitialized(false)
1876 , shapes(other.shapes)
1877 , shapes_faces(other.shapes_faces)
1878 , is_interior(other.is_interior)
1879{
1880 connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1881 [this]() { this->is_reinitialized = false; });
1882}
1883
1884
1885
1886template <int n_components_, int dim, int spacedim, typename Number>
1889{
1890 connection_is_reinitialized.disconnect();
1891}
1892
1893
1894
1895template <int n_components_, int dim, int spacedim, typename Number>
1896void
1898 const unsigned int first_selected_component)
1899{
1900 AssertIndexRange(first_selected_component + n_components,
1901 fe->n_components() + 1);
1902
1903 shapes.reserve(100);
1904
1905 bool same_base_element = true;
1906 unsigned int base_element_number = 0;
1907 component_in_base_element = 0;
1908 unsigned int component = 0;
1909 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1910 if (component + fe->element_multiplicity(base_element_number) >
1911 first_selected_component)
1912 {
1913 if (first_selected_component + n_components >
1914 component + fe->element_multiplicity(base_element_number))
1915 same_base_element = false;
1916 component_in_base_element = first_selected_component - component;
1917 break;
1918 }
1919 else
1920 component += fe->element_multiplicity(base_element_number);
1921
1924 *fe, base_element_number) &&
1925 same_base_element)
1926 {
1927 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1928 renumber = shape_info.lexicographic_numbering;
1929 dofs_per_component = shape_info.dofs_per_component_on_cell;
1930 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1932 fe->base_element(base_element_number));
1933
1934 bool is_lexicographic = true;
1935 for (unsigned int i = 0; i < renumber.size(); ++i)
1936 if (i != renumber[i])
1937 is_lexicographic = false;
1938
1939 if (is_lexicographic)
1940 renumber.clear();
1941
1942 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1943 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1944 poly[1].value(1.) == 1.) &&
1945 (fe->n_components() == n_components);
1946
1947 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1948 const unsigned int size_cell = dofs_per_component * n_components;
1949 scratch_data_scalar.resize(size_face + size_cell);
1950
1951 solution_renumbered.resize(dofs_per_component);
1952 solution_renumbered_vectorized.resize(dofs_per_component);
1953
1954 fast_path = true;
1955 }
1956 else
1957 {
1958 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1959 for (unsigned int d = 0; d < n_components; ++d)
1960 {
1961 const unsigned int component = first_selected_component + d;
1962 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1963 {
1964 const bool is_primitive =
1965 fe->is_primitive() || fe->is_primitive(i);
1966 if (is_primitive)
1967 nonzero_shape_function_component[i][d] =
1968 (component == fe->system_to_component_index(i).first);
1969 else
1970 nonzero_shape_function_component[i][d] =
1971 (fe->get_nonzero_components(i)[component] == true);
1972 }
1973 }
1974
1975 fast_path = false;
1976 }
1977}
1978
1979
1980
1981template <int n_components_, int dim, int spacedim, typename Number>
1982template <bool is_face, bool is_linear>
1983inline void
1985{
1986 const unsigned int geometry_index =
1987 mapping_info->template compute_geometry_index_offset<is_face>(
1988 current_cell_index, current_face_number);
1989
1990 cell_type = mapping_info->get_cell_type(geometry_index);
1991
1992 const_cast<unsigned int &>(n_q_points_scalar) =
1993 mapping_info->get_n_q_points_unvectorized(geometry_index);
1994
1995 // round up n_q_points_scalar / n_lanes_internal
1996 const_cast<unsigned int &>(n_q_batches) =
1997 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
1998
1999 const unsigned int n_q_points_before = n_q_points;
2000
2001 const_cast<unsigned int &>(n_q_points) =
2002 (stride == 1) ? n_q_batches : n_q_points_scalar;
2003
2004 if (n_q_points != n_q_points_before)
2005 {
2006 if (update_flags & update_values)
2007 values.resize(n_q_points);
2008 if (update_flags & update_gradients)
2009 gradients.resize(n_q_points);
2010 }
2011
2012 if (n_q_points == 0)
2013 {
2014 is_reinitialized = true;
2015 return;
2016 }
2017
2018 // set unit point pointer
2019 const unsigned int unit_point_offset =
2020 mapping_info->compute_unit_point_index_offset(geometry_index);
2021
2022 if (is_face)
2023 unit_point_faces_ptr =
2024 mapping_info->get_unit_point_faces(unit_point_offset);
2025 else
2026 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2027
2028 // set data pointers
2029 const unsigned int data_offset =
2030 mapping_info->compute_data_index_offset(geometry_index);
2031 const unsigned int compressed_data_offset =
2032 mapping_info->compute_compressed_data_index_offset(geometry_index);
2033#ifdef DEBUG
2034 const UpdateFlags update_flags_mapping =
2035 mapping_info->get_update_flags_mapping();
2036 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2037 real_point_ptr = mapping_info->get_real_point(data_offset);
2038 if (update_flags_mapping & UpdateFlags::update_jacobians)
2039 jacobian_ptr =
2040 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2041 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2042 inverse_jacobian_ptr =
2043 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2044 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2045 normal_ptr = mapping_info->get_normal_vector(data_offset);
2046 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2047 JxW_ptr = mapping_info->get_JxW(data_offset);
2048#else
2049 real_point_ptr = mapping_info->get_real_point(data_offset);
2050 jacobian_ptr =
2051 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2052 inverse_jacobian_ptr =
2053 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2054 normal_ptr = mapping_info->get_normal_vector(data_offset);
2055 JxW_ptr = mapping_info->get_JxW(data_offset);
2056#endif
2057
2058 if (!is_linear && fast_path)
2059 {
2060 const std::size_t n_shapes = poly.size();
2061
2062 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2063 if (is_face)
2064 {
2065 if (dim > 1)
2066 {
2067 shapes_faces.resize_fast(n_q_batches * n_shapes);
2069 shapes_faces.data() + qb * n_shapes,
2070 poly,
2071 unit_point_faces_ptr[qb],
2072 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2073 }
2074 }
2075 else
2076 {
2077 shapes.resize_fast(n_q_batches * n_shapes);
2079 shapes.data() + qb * n_shapes,
2080 poly,
2081 unit_point_ptr[qb],
2082 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2083 }
2084 }
2085
2086 is_reinitialized = true;
2087}
2088
2089
2090
2091template <int n_components_, int dim, int spacedim, typename Number>
2092inline const typename FEPointEvaluationBase<n_components_,
2093 dim,
2094 spacedim,
2095 Number>::value_type &
2097 const unsigned int point_index) const
2098{
2099 AssertIndexRange(point_index, values.size());
2100 return values[point_index];
2101}
2102
2103
2104
2105template <int n_components_, int dim, int spacedim, typename Number>
2106inline const typename FEPointEvaluationBase<n_components_,
2107 dim,
2108 spacedim,
2109 Number>::gradient_type &
2111 const unsigned int point_index) const
2112{
2113 AssertIndexRange(point_index, gradients.size());
2114 return gradients[point_index];
2115}
2116
2117
2118
2119template <int n_components_, int dim, int spacedim, typename Number>
2120inline void
2122 const value_type &value,
2123 const unsigned int point_index)
2124{
2125 AssertIndexRange(point_index, n_q_points);
2126 values[point_index] = value;
2127}
2128
2129
2130
2131template <int n_components_, int dim, int spacedim, typename Number>
2132inline void
2134 const gradient_type &gradient,
2135 const unsigned int point_index)
2136{
2137 AssertIndexRange(point_index, n_q_points);
2138 gradients[point_index] = gradient;
2139}
2140
2141
2142
2143template <int n_components_, int dim, int spacedim, typename Number>
2146 const unsigned int point_index) const
2147{
2148 AssertIndexRange(point_index, n_q_points);
2149 Assert(jacobian_ptr != nullptr,
2151 ExcFEPointEvaluationAccessToUninitializedMappingField(
2152 "update_jacobians"));
2153 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2154 GeometryType::affine ?
2155 0 :
2156 point_index];
2157}
2158
2159
2160
2161template <int n_components_, int dim, int spacedim, typename Number>
2164 const unsigned int point_index) const
2165{
2166 AssertIndexRange(point_index, n_q_points);
2167 Assert(inverse_jacobian_ptr != nullptr,
2169 ExcFEPointEvaluationAccessToUninitializedMappingField(
2170 "update_inverse_jacobians"));
2171 return inverse_jacobian_ptr
2172 [cell_type <=
2174 0 :
2175 point_index];
2176}
2177
2178
2179
2180template <int n_components_, int dim, int spacedim, typename Number>
2181inline Number
2183 const unsigned int point_index) const
2184{
2185 AssertIndexRange(point_index, n_q_points);
2186 Assert(JxW_ptr != nullptr,
2188 ExcFEPointEvaluationAccessToUninitializedMappingField(
2189 "update_JxW_values"));
2190 return JxW_ptr[point_index];
2191}
2192
2193
2194
2195template <int n_components_, int dim, int spacedim, typename Number>
2198 const unsigned int point_index) const
2199{
2200 return quadrature_point(point_index);
2201 ;
2202}
2203
2204
2205
2206template <int n_components_, int dim, int spacedim, typename Number>
2209 const unsigned int point_index) const
2210{
2211 AssertIndexRange(point_index, n_q_points);
2212 Assert(real_point_ptr != nullptr,
2214 ExcFEPointEvaluationAccessToUninitializedMappingField(
2215 "update_quadrature_points"));
2216 return real_point_ptr[point_index];
2217}
2218
2219
2220
2221template <int n_components_, int dim, int spacedim, typename Number>
2222inline Point<dim, Number>
2224 const unsigned int point_index) const
2225{
2226 AssertIndexRange(point_index, n_q_points);
2227 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2228 Point<dim, Number> unit_point;
2229 for (unsigned int d = 0; d < dim; ++d)
2231 unit_point_ptr[point_index / stride][d], point_index % stride);
2232 return unit_point;
2233}
2234
2235
2236
2237template <int n_components_, int dim, int spacedim, typename Number>
2244
2245
2246
2247template <int n_components_, int dim, int spacedim, typename Number>
2251 const unsigned int first_selected_component)
2252 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2253 mapping_info,
2254 fe,
2255 first_selected_component)
2256{}
2257
2258
2259
2260template <int n_components_, int dim, int spacedim, typename Number>
2262 const Mapping<dim, spacedim> &mapping,
2264 const UpdateFlags update_flags,
2265 const unsigned int first_selected_component)
2266 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2267 mapping,
2268 fe,
2269 update_flags,
2270 first_selected_component)
2271{}
2272
2273
2274
2275template <int n_components_, int dim, int spacedim, typename Number>
2276inline void
2278{
2279 this->current_cell_index = numbers::invalid_unsigned_int;
2280 this->current_face_number = numbers::invalid_unsigned_int;
2281
2282 if (this->use_linear_path)
2283 this->template do_reinit<false, true>();
2284 else
2285 this->template do_reinit<false, false>();
2286}
2287
2288
2289
2290template <int n_components_, int dim, int spacedim, typename Number>
2291inline void
2294 const ArrayView<const Point<dim>> &unit_points)
2295{
2296 // reinit is only allowed for mapping computation on the fly
2297 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2299
2300 this->mapping_info->reinit(cell, unit_points);
2301
2302 if (!this->fast_path)
2303 {
2304 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2305 *this->mapping,
2306 *this->fe,
2308 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2309 this->update_flags);
2310 this->fe_values->reinit(cell);
2311 }
2312
2313 if (this->use_linear_path)
2314 this->template do_reinit<false, true>();
2315 else
2316 this->template do_reinit<false, false>();
2317}
2318
2319
2320
2321template <int n_components_, int dim, int spacedim, typename Number>
2322inline void
2324 const unsigned int cell_index)
2325{
2326 this->current_cell_index = cell_index;
2327 this->current_face_number = numbers::invalid_unsigned_int;
2328
2329 if (this->use_linear_path)
2330 this->template do_reinit<false, true>();
2331 else
2332 this->template do_reinit<false, false>();
2333
2334 if (!this->fast_path)
2335 {
2336 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2337
2338 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2339 for (unsigned int d = 0; d < dim; ++d)
2340 unit_points[v][d] =
2341 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2342
2343 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2344 *this->mapping,
2345 *this->fe,
2347 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2348 this->update_flags);
2349
2350 this->fe_values->reinit(
2351 this->mapping_info->get_cell_iterator(this->current_cell_index));
2352 }
2353}
2354
2355
2356
2357template <int n_components_, int dim, int spacedim, typename Number>
2358template <std::size_t stride_view>
2359void
2362 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2363{
2364 if (!this->is_reinitialized)
2365 reinit();
2366
2367 if (this->n_q_points == 0)
2368 return;
2369
2370 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2371
2372 if (!((evaluation_flags & EvaluationFlags::values) ||
2373 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2374 return;
2375
2376 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2377 if (this->fast_path)
2378 {
2379 if (this->use_linear_path)
2380 evaluate_fast<true>(solution_values, evaluation_flags);
2381 else
2382 evaluate_fast<false>(solution_values, evaluation_flags);
2383 }
2384 else
2385 evaluate_slow(solution_values, evaluation_flags);
2386}
2387
2388
2389
2390template <int n_components_, int dim, int spacedim, typename Number>
2391void
2393 const ArrayView<const ScalarNumber> &solution_values,
2394 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2395{
2396 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2397 solution_values.size()),
2398 evaluation_flags);
2399}
2400
2401
2402
2403template <int n_components_, int dim, int spacedim, typename Number>
2404template <std::size_t stride_view>
2405void
2407 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2408 const EvaluationFlags::EvaluationFlags &integration_flags,
2409 const bool sum_into_values)
2410{
2411 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2412}
2413
2414
2415
2416template <int n_components_, int dim, int spacedim, typename Number>
2417void
2419 const ArrayView<ScalarNumber> &solution_values,
2420 const EvaluationFlags::EvaluationFlags &integration_flags,
2421 const bool sum_into_values)
2422{
2423 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2424 solution_values.size()),
2425 integration_flags,
2426 sum_into_values);
2427}
2428
2429
2430
2431template <int n_components_, int dim, int spacedim, typename Number>
2433 scalar_value_type
2435 const
2436{
2437 value_type return_value = {};
2438
2439 for (const auto point_index : this->quadrature_point_indices())
2440 return_value += values[point_index] * this->JxW(point_index);
2441
2442 return ETT::sum_value(return_value);
2443}
2444
2445
2446
2447template <int n_components_, int dim, int spacedim, typename Number>
2448template <std::size_t stride_view>
2449void
2451 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2452 const EvaluationFlags::EvaluationFlags &integration_flags,
2453 const bool sum_into_values)
2454{
2455 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2456}
2457
2458
2459
2460template <int n_components_, int dim, int spacedim, typename Number>
2461void
2463 const ArrayView<ScalarNumber> &solution_values,
2464 const EvaluationFlags::EvaluationFlags &integration_flags,
2465 const bool sum_into_values)
2466{
2467 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2468 solution_values.size()),
2469 integration_flags,
2470 sum_into_values);
2471}
2472
2473
2474
2475template <int n_components_, int dim, int spacedim, typename Number>
2476template <bool is_linear, std::size_t stride_view>
2477inline void
2480{
2481 const unsigned int dofs_per_comp =
2482 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2483
2484 for (unsigned int comp = 0; comp < n_components; ++comp)
2485 {
2486 const std::size_t offset =
2487 (this->component_in_base_element + comp) * dofs_per_comp;
2488
2489 if ((is_linear && n_components == 1) || this->renumber.empty())
2490 {
2491 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2492 ETT::read_value(solution_values[i + offset],
2493 comp,
2494 this->solution_renumbered[i]);
2495 }
2496 else
2497 {
2498 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2499 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2500 ETT::read_value(solution_values[renumber_ptr[i]],
2501 comp,
2502 this->solution_renumbered[i]);
2503 }
2504 }
2505}
2506
2507
2508
2509template <int n_components_, int dim, int spacedim, typename Number>
2510template <bool is_linear, std::size_t stride_view>
2511inline void
2514 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2515 const unsigned int n_shapes,
2516 const unsigned int qb,
2517 vectorized_value_type &value,
2519{
2520 if (evaluation_flags & EvaluationFlags::gradients)
2521 {
2522 std::array<vectorized_value_type, dim + 1> result;
2523 if constexpr (is_linear)
2524 {
2525 if constexpr (n_components == 1)
2526 result =
2528 dim,
2531 1,
2532 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2533 else
2534 result =
2536 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2537 }
2538 else
2540 dim,
2543 1,
2544 false>(this->shapes.data() + qb * n_shapes,
2545 n_shapes,
2546 this->solution_renumbered.data());
2547 gradient[0] = result[0];
2548 if (dim > 1)
2549 gradient[1] = result[1];
2550 if (dim > 2)
2551 gradient[2] = result[2];
2552 value = result[dim];
2553 }
2554 else
2555 {
2556 if constexpr (is_linear)
2557 {
2558 if constexpr (n_components == 1)
2560 dim,
2563 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2564 else
2566 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2567 }
2568 else
2569 value =
2573 false>(
2574 this->shapes.data() + qb * n_shapes,
2575 n_shapes,
2576 this->solution_renumbered.data());
2577 }
2578}
2579
2580
2581
2582template <int n_components_, int dim, int spacedim, typename Number>
2583template <bool is_linear, std::size_t stride_view>
2584inline void
2587 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2588{
2589 if (!(is_linear && n_components == 1))
2590 prepare_evaluate_fast<is_linear>(solution_values);
2591
2592 // loop over quadrature batches qb
2593 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2594
2595 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2596 {
2599
2600 compute_evaluate_fast<is_linear>(
2601 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2602
2603 if (evaluation_flags & EvaluationFlags::values)
2604 {
2605 for (unsigned int v = 0, offset = qb * stride;
2606 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2607 ++v, ++offset)
2608 ETT::set_value(value, v, this->values[offset]);
2609 }
2610 if (evaluation_flags & EvaluationFlags::gradients)
2611 {
2612 Assert(this->update_flags & update_gradients ||
2613 this->update_flags & update_inverse_jacobians,
2615
2616 for (unsigned int v = 0, offset = qb * stride;
2617 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2618 ++v, ++offset)
2619 {
2620 unit_gradient_type unit_gradient;
2621 ETT::set_gradient(gradient, v, unit_gradient);
2622 this->gradients[offset] =
2623 this->cell_type <=
2626 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2628 this
2629 ->inverse_jacobian_ptr[this->cell_type <=
2631 GeometryType::affine ?
2632 0 :
2633 offset]
2634 .transpose(),
2635 unit_gradient);
2636 }
2637 }
2638 }
2639}
2640
2641
2642
2643template <int n_components_, int dim, int spacedim, typename Number>
2644template <std::size_t stride_view>
2645inline void
2648 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2649{
2650 // slow path with FEValues
2651 Assert(this->fe_values.get() != nullptr,
2652 ExcMessage(
2653 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2654
2655 const std::size_t n_points = this->fe_values->get_quadrature().size();
2656
2657 if (evaluation_flags & EvaluationFlags::values)
2658 {
2659 this->values.resize(this->n_q_points);
2660 std::fill(this->values.begin(), this->values.end(), value_type());
2661 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2662 {
2663 const ScalarNumber value = solution_values[i];
2664 for (unsigned int d = 0; d < n_components; ++d)
2665 if (this->nonzero_shape_function_component[i][d] &&
2666 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2667 for (unsigned int qb = 0, q = 0; q < n_points;
2668 ++qb, q += n_lanes_user_interface)
2669 for (unsigned int v = 0;
2670 v < n_lanes_user_interface && q + v < n_points;
2671 ++v)
2672 ETT::access(this->values[qb],
2673 v,
2674 d,
2675 this->fe_values->shape_value(i, q + v) * value);
2676 else if (this->nonzero_shape_function_component[i][d])
2677 for (unsigned int qb = 0, q = 0; q < n_points;
2678 ++qb, q += n_lanes_user_interface)
2679 for (unsigned int v = 0;
2680 v < n_lanes_user_interface && q + v < n_points;
2681 ++v)
2682 ETT::access(this->values[qb],
2683 v,
2684 d,
2685 this->fe_values->shape_value_component(i,
2686 q + v,
2687 d) *
2688 value);
2689 }
2690 }
2691
2692 if (evaluation_flags & EvaluationFlags::gradients)
2693 {
2694 this->gradients.resize(this->n_q_points);
2695 std::fill(this->gradients.begin(),
2696 this->gradients.end(),
2697 gradient_type());
2698 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2699 {
2700 const ScalarNumber value = solution_values[i];
2701 for (unsigned int d = 0; d < n_components; ++d)
2702 if (this->nonzero_shape_function_component[i][d] &&
2703 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2704 for (unsigned int qb = 0, q = 0; q < n_points;
2705 ++qb, q += n_lanes_user_interface)
2706 for (unsigned int v = 0;
2707 v < n_lanes_user_interface && q + v < n_points;
2708 ++v)
2709 ETT::access(this->gradients[qb],
2710 v,
2711 d,
2712 this->fe_values->shape_grad(i, q + v) * value);
2713 else if (this->nonzero_shape_function_component[i][d])
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(
2720 this->gradients[qb],
2721 v,
2722 d,
2723 this->fe_values->shape_grad_component(i, q + v, d) * value);
2724 }
2725 }
2726}
2727
2728
2729
2730template <int n_components_, int dim, int spacedim, typename Number>
2731template <bool is_linear>
2732inline void
2734 const EvaluationFlags::EvaluationFlags &integration_flags,
2735 const unsigned int n_shapes,
2736 const unsigned int qb,
2737 const vectorized_value_type value,
2739 vectorized_value_type *solution_values_vectorized_linear)
2740{
2741 if (integration_flags & EvaluationFlags::gradients)
2743 is_linear,
2744 dim,
2746 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2747 n_shapes,
2748 &value,
2749 gradient,
2750 is_linear ?
2751 solution_values_vectorized_linear :
2752 this->solution_renumbered_vectorized.data(),
2753 this->unit_point_ptr[qb],
2754 qb != 0);
2755 else
2757 dim,
2760 this->shapes.data() + qb * n_shapes,
2761 n_shapes,
2762 value,
2763 is_linear ? solution_values_vectorized_linear :
2764 this->solution_renumbered_vectorized.data(),
2765 this->unit_point_ptr[qb],
2766 qb != 0);
2767}
2768
2769
2770
2771template <int n_components_, int dim, int spacedim, typename Number>
2772template <bool is_linear, std::size_t stride_view>
2773inline void
2775 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2776 vectorized_value_type *solution_values_vectorized_linear,
2777 const bool sum_into_values)
2778{
2779 if (!sum_into_values && this->fe->n_components() > n_components)
2780 for (unsigned int i = 0; i < solution_values.size(); ++i)
2781 solution_values[i] = 0;
2782
2783 const unsigned int dofs_per_comp =
2784 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2785
2786 for (unsigned int comp = 0; comp < n_components; ++comp)
2787 {
2788 const std::size_t offset =
2789 (this->component_in_base_element + comp) * dofs_per_comp;
2790
2791 if (is_linear || this->renumber.empty())
2792 {
2793 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2794 if (sum_into_values)
2795 solution_values[i + offset] +=
2796 ETT::sum_value(comp,
2797 is_linear ?
2798 *(solution_values_vectorized_linear + i) :
2799 this->solution_renumbered_vectorized[i]);
2800 else
2801 solution_values[i + offset] =
2802 ETT::sum_value(comp,
2803 is_linear ?
2804 *(solution_values_vectorized_linear + i) :
2805 this->solution_renumbered_vectorized[i]);
2806 }
2807 else
2808 {
2809 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2810 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2811 if (sum_into_values)
2812 solution_values[renumber_ptr[i]] +=
2813 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2814 else
2815 solution_values[renumber_ptr[i]] =
2816 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2817 }
2818 }
2819}
2820
2821
2822
2823template <int n_components_, int dim, int spacedim, typename Number>
2824template <bool do_JxW, bool is_linear, std::size_t stride_view>
2825inline void
2827 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2828 const EvaluationFlags::EvaluationFlags &integration_flags,
2829 const bool sum_into_values)
2830{
2831 // zero out lanes of incomplete last quadrature point batch
2832 if constexpr (stride == 1)
2833 if (const unsigned int n_filled_lanes =
2834 this->n_q_points_scalar & (n_lanes_internal - 1);
2835 n_filled_lanes > 0)
2836 {
2837 if (integration_flags & EvaluationFlags::values)
2838 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2839 ETT::set_zero_value(this->values.back(), v);
2840 if (integration_flags & EvaluationFlags::gradients)
2841 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2842 ETT::set_zero_gradient(this->gradients.back(), v);
2843 }
2844
2845 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
2846 solution_values_vectorized_linear = {};
2847
2848 // loop over quadrature batches qb
2849 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2850
2851 const bool cartesian_cell =
2853 const bool affine_cell =
2855 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2856 {
2857 vectorized_value_type value = {};
2859
2860 if (integration_flags & EvaluationFlags::values)
2861 for (unsigned int v = 0, offset = qb * stride;
2862 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2863 ++v, ++offset)
2864 ETT::get_value(value,
2865 v,
2866 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
2867 this->values[offset]);
2868
2869 if (integration_flags & EvaluationFlags::gradients)
2870 for (unsigned int v = 0, offset = qb * stride;
2871 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2872 ++v, ++offset)
2873 {
2874 const gradient_type grad_w =
2875 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
2876 this->gradients[offset];
2877 ETT::get_gradient(
2878 gradient,
2879 v,
2880 cartesian_cell ?
2881 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
2882 grad_w) :
2884 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
2885 grad_w));
2886 }
2887
2888 compute_integrate_fast<is_linear>(
2889 integration_flags,
2890 n_shapes,
2891 qb,
2892 value,
2893 gradient,
2894 solution_values_vectorized_linear.data());
2895 }
2896
2897 // add between the lanes and write into the result
2898 finish_integrate_fast<is_linear>(solution_values,
2899 solution_values_vectorized_linear.data(),
2900 sum_into_values);
2901}
2902
2903
2904
2905template <int n_components_, int dim, int spacedim, typename Number>
2906template <bool do_JxW, std::size_t stride_view>
2907inline void
2909 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2910 const EvaluationFlags::EvaluationFlags &integration_flags,
2911 const bool sum_into_values)
2912{
2913 // slow path with FEValues
2914 Assert(this->fe_values.get() != nullptr,
2915 ExcMessage(
2916 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2917 if (!sum_into_values)
2918 for (unsigned int i = 0; i < solution_values.size(); ++i)
2919 solution_values[i] = 0;
2920
2921 const std::size_t n_points = this->fe_values->get_quadrature().size();
2922
2923 if (integration_flags & EvaluationFlags::values)
2924 {
2925 AssertIndexRange(this->n_q_points, this->values.size() + 1);
2926 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2927 {
2928 for (unsigned int d = 0; d < n_components; ++d)
2929 if (this->nonzero_shape_function_component[i][d] &&
2930 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2931 for (unsigned int qb = 0, q = 0; q < n_points;
2932 ++qb, q += n_lanes_user_interface)
2933 for (unsigned int v = 0;
2934 v < n_lanes_user_interface && q + v < n_points;
2935 ++v)
2936 solution_values[i] +=
2937 this->fe_values->shape_value(i, q + v) *
2938 ETT::access(this->values[qb], v, d) *
2939 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2940 else if (this->nonzero_shape_function_component[i][d])
2941 for (unsigned int qb = 0, q = 0; q < n_points;
2942 ++qb, q += n_lanes_user_interface)
2943 for (unsigned int v = 0;
2944 v < n_lanes_user_interface && q + v < n_points;
2945 ++v)
2946 solution_values[i] +=
2947 this->fe_values->shape_value_component(i, q + v, d) *
2948 ETT::access(this->values[qb], v, d) *
2949 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2950 }
2951 }
2952
2953 if (integration_flags & EvaluationFlags::gradients)
2954 {
2955 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
2956 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2957 {
2958 for (unsigned int d = 0; d < n_components; ++d)
2959 if (this->nonzero_shape_function_component[i][d] &&
2960 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2961 for (unsigned int qb = 0, q = 0; q < n_points;
2962 ++qb, q += n_lanes_user_interface)
2963 for (unsigned int v = 0;
2964 v < n_lanes_user_interface && q + v < n_points;
2965 ++v)
2966 solution_values[i] +=
2967 this->fe_values->shape_grad(i, q + v) *
2968 ETT::access(this->gradients[qb], v, d) *
2969 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2970 else if (this->nonzero_shape_function_component[i][d])
2971 for (unsigned int qb = 0, q = 0; q < n_points;
2972 ++qb, q += n_lanes_user_interface)
2973 for (unsigned int v = 0;
2974 v < n_lanes_user_interface && q + v < n_points;
2975 ++v)
2976 solution_values[i] +=
2977 this->fe_values->shape_grad_component(i, q + v, d) *
2978 ETT::access(this->gradients[qb], v, d) *
2979 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2980 }
2981 }
2982}
2983
2984
2985
2986template <int n_components_, int dim, int spacedim, typename Number>
2987template <bool do_JxW, std::size_t stride_view>
2988void
2990 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2991 const EvaluationFlags::EvaluationFlags &integration_flags,
2992 const bool sum_into_values)
2993{
2994 if (!this->is_reinitialized)
2995 reinit();
2996
2997 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
2998
2999 if (this->n_q_points == 0 || // no evaluation points provided
3000 !((integration_flags & EvaluationFlags::values) ||
3001 (integration_flags &
3002 EvaluationFlags::gradients))) // no integration flags
3003 {
3004 if (!sum_into_values)
3005 for (unsigned int i = 0; i < solution_values.size(); ++i)
3006 solution_values[i] = 0;
3007 return;
3008 }
3009
3010 Assert(
3011 !do_JxW || this->JxW_ptr != nullptr,
3012 ExcMessage(
3013 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3014
3015 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3016 if (this->fast_path)
3017 {
3018 if (this->use_linear_path)
3019 integrate_fast<do_JxW, true>(solution_values,
3020 integration_flags,
3021 sum_into_values);
3022 else
3023 integrate_fast<do_JxW, false>(solution_values,
3024 integration_flags,
3025 sum_into_values);
3026 }
3027 else
3028 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3029}
3030
3031
3032
3033template <int n_components_, int dim, int spacedim, typename Number>
3036 const unsigned int point_index) const
3037{
3038 AssertIndexRange(point_index, this->n_q_points);
3039 Assert(this->normal_ptr != nullptr,
3041 ExcFEPointEvaluationAccessToUninitializedMappingField(
3042 "update_normal_vectors"));
3043 if (this->is_interior)
3044 return this->normal_ptr[point_index];
3045 else
3046 return -this->normal_ptr[point_index];
3047}
3048
3049
3050
3051template <int n_components_, int dim, int spacedim, typename Number>
3056 const bool is_interior,
3057 const unsigned int first_selected_component)
3058 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3059 mapping_info,
3060 fe,
3061 first_selected_component,
3062 is_interior)
3063{}
3064
3065
3066
3067template <int n_components_, int dim, int spacedim, typename Number>
3068inline void
3070 const unsigned int cell_index,
3071 const unsigned int face_number)
3072{
3073 this->current_cell_index = cell_index;
3074 this->current_face_number = face_number;
3075
3076 if (this->use_linear_path)
3077 this->template do_reinit<true, true>();
3078 else
3079 this->template do_reinit<true, false>();
3080}
3081
3082
3083
3084template <int n_components_, int dim, int spacedim, typename Number>
3085inline void
3087 const unsigned int face_index)
3088{
3089 this->current_cell_index = face_index;
3090 this->current_face_number =
3091 this->mapping_info->get_face_number(face_index, this->is_interior);
3092
3093 if (this->use_linear_path)
3094 this->template do_reinit<true, true>();
3095 else
3096 this->template do_reinit<true, false>();
3097}
3098
3099
3100
3101template <int n_components_, int dim, int spacedim, typename Number>
3102template <std::size_t stride_view>
3103void
3106 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3107{
3108 Assert(this->is_reinitialized, ExcMessage("Is not reinitialized!"));
3109
3110 if (this->n_q_points == 0)
3111 return;
3112
3113 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3114
3115 if (!((evaluation_flags & EvaluationFlags::values) ||
3116 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3117 return;
3118
3119 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3120
3121 if (this->use_linear_path)
3122 do_evaluate<true>(solution_values, evaluation_flags);
3123 else
3124 do_evaluate<false>(solution_values, evaluation_flags);
3125}
3126
3127
3128
3129template <int n_components_, int dim, int spacedim, typename Number>
3130void
3132 const ArrayView<const ScalarNumber> &solution_values,
3133 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3134{
3135 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3136 solution_values.size()),
3137 evaluation_flags);
3138}
3139
3140
3141
3142template <int n_components_, int dim, int spacedim, typename Number>
3143template <bool is_linear, std::size_t stride_view>
3144void
3147 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3148{
3149 const unsigned int dofs_per_comp =
3150 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3151
3152 const ScalarNumber *input;
3153 if (stride_view == 1 && this->component_in_base_element == 0 &&
3154 (is_linear || this->renumber.empty()))
3155 input = solution_values.data();
3156 else
3157 {
3158 for (unsigned int comp = 0; comp < n_components; ++comp)
3159 {
3160 const std::size_t offset =
3161 (this->component_in_base_element + comp) * dofs_per_comp;
3162
3163 if (is_linear || this->renumber.empty())
3164 {
3165 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3166 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3167 solution_values[i + offset];
3168 }
3169 else
3170 {
3171 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3172 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3173 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3174 solution_values[renumber_ptr[i]];
3175 }
3176 }
3177 input = this->scratch_data_scalar.data();
3178 }
3179
3180 ScalarNumber *output =
3181 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3182
3183 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3184 template interpolate<true, false>(n_components,
3185 evaluation_flags,
3186 this->shape_info,
3187 input,
3188 output,
3189 this->current_face_number);
3190
3191 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3192}
3193
3194
3195
3196template <int n_components_, int dim, int spacedim, typename Number>
3197template <std::size_t stride_view>
3198void
3200 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3201 const EvaluationFlags::EvaluationFlags &integration_flags,
3202 const bool sum_into_values)
3203{
3204 Assert(this->is_reinitialized, ExcMessage("Is not reinitialized!"));
3205
3206 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3207
3208 if (this->n_q_points == 0 || // no evaluation points provided
3209 !((integration_flags & EvaluationFlags::values) ||
3210 (integration_flags &
3211 EvaluationFlags::gradients))) // no integration flags
3212 {
3213 if (!sum_into_values)
3214 for (unsigned int i = 0; i < solution_values.size(); ++i)
3215 solution_values[i] = 0;
3216 return;
3217 }
3218
3219 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3220
3221 if (this->use_linear_path)
3222 do_integrate<true, true>(solution_values,
3223 integration_flags,
3224 sum_into_values);
3225 else
3226 do_integrate<true, false>(solution_values,
3227 integration_flags,
3228 sum_into_values);
3229}
3230
3231
3232
3233template <int n_components_, int dim, int spacedim, typename Number>
3234void
3236 const ArrayView<ScalarNumber> &solution_values,
3237 const EvaluationFlags::EvaluationFlags &integration_flags,
3238 const bool sum_into_values)
3239{
3240 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3241 solution_values.size()),
3242 integration_flags,
3243 sum_into_values);
3244}
3245
3246
3247
3248template <int n_components_, int dim, int spacedim, typename Number>
3249template <std::size_t stride_view>
3250void
3252 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3253 const EvaluationFlags::EvaluationFlags &integration_flags,
3254 const bool sum_into_values)
3255{
3256 Assert(this->is_reinitialized, ExcMessage("Is not reinitialized!"));
3257
3258 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3259
3260 if (this->n_q_points == 0 || // no evaluation points provided
3261 !((integration_flags & EvaluationFlags::values) ||
3262 (integration_flags &
3263 EvaluationFlags::gradients))) // no integration flags
3264 {
3265 if (!sum_into_values)
3266 for (unsigned int i = 0; i < solution_values.size(); ++i)
3267 solution_values[i] = 0;
3268 return;
3269 }
3270
3271 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3272
3273 if (this->use_linear_path)
3274 do_integrate<false, true>(solution_values,
3275 integration_flags,
3276 sum_into_values);
3277 else
3278 do_integrate<false, false>(solution_values,
3279 integration_flags,
3280 sum_into_values);
3281}
3282
3283
3284
3285template <int n_components_, int dim, int spacedim, typename Number>
3286void
3288 const ArrayView<ScalarNumber> &solution_values,
3289 const EvaluationFlags::EvaluationFlags &integration_flags,
3290 const bool sum_into_values)
3291{
3292 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3293 solution_values.size()),
3294 integration_flags,
3295 sum_into_values);
3296}
3297
3298
3299
3300template <int n_components_, int dim, int spacedim, typename Number>
3301template <bool do_JxW, bool is_linear, std::size_t stride_view>
3302void
3304 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3305 const EvaluationFlags::EvaluationFlags &integration_flags,
3306 const bool sum_into_values)
3307{
3308 if (!sum_into_values && this->fe->n_components() > n_components)
3309 for (unsigned int i = 0; i < solution_values.size(); ++i)
3310 solution_values[i] = 0;
3311
3312 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3313 integration_flags,
3314 false);
3315
3316 ScalarNumber *input = this->scratch_data_scalar.begin();
3317
3318 if (stride_view == 1 && this->component_in_base_element == 0 &&
3319 (is_linear || this->renumber.empty()))
3320 {
3321 if (sum_into_values)
3322 internal::
3323 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3324 template interpolate<false, true>(n_components,
3325 integration_flags,
3326 this->shape_info,
3327 input,
3328 solution_values.data(),
3329 this->current_face_number);
3330 else
3331 internal::
3332 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3333 template interpolate<false, false>(n_components,
3334 integration_flags,
3335 this->shape_info,
3336 input,
3337 solution_values.data(),
3338 this->current_face_number);
3339 }
3340 else
3341 {
3342 const unsigned int dofs_per_comp_face =
3343 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3344
3345 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3346 ScalarNumber *output = input + size_input;
3347
3348 internal::
3349 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3350 template interpolate<false, false>(n_components,
3351 integration_flags,
3352 this->shape_info,
3353 input,
3354 output,
3355 this->current_face_number);
3356
3357 const unsigned int dofs_per_comp =
3358 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3359
3360 for (unsigned int comp = 0; comp < n_components; ++comp)
3361 {
3362 const std::size_t offset =
3363 (this->component_in_base_element + comp) * dofs_per_comp;
3364
3365 if (is_linear || this->renumber.empty())
3366 {
3367 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3368 if (sum_into_values)
3369 solution_values[i + offset] +=
3370 output[i + comp * dofs_per_comp];
3371 else
3372 solution_values[i + offset] =
3373 output[i + comp * dofs_per_comp];
3374 }
3375 else
3376 {
3377 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3378 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3379 if (sum_into_values)
3380 solution_values[renumber_ptr[i]] +=
3381 output[i + comp * dofs_per_comp];
3382 else
3383 solution_values[renumber_ptr[i]] =
3384 output[i + comp * dofs_per_comp];
3385 }
3386 }
3387 }
3388}
3389
3390
3391
3392template <int n_components_, int dim, int spacedim, typename Number>
3393template <int stride_face_dof>
3394void
3396 const ScalarNumber *face_dof_values,
3397 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3398{
3399 if (this->use_linear_path)
3400 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3401 evaluation_flags);
3402 else
3403 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3404 evaluation_flags);
3405}
3406
3407
3408
3409template <int n_components_, int dim, int spacedim, typename Number>
3410template <bool is_linear, int stride_face_dof>
3411inline void
3413 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3414 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3415{
3416 const scalar_value_type *face_dof_values_ptr;
3417 if constexpr (n_components == 1)
3418 face_dof_values_ptr = face_dof_values;
3419 else
3420 {
3421 const unsigned int dofs_per_comp_face =
3422 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3423 for (unsigned int comp = 0; comp < n_components; ++comp)
3424 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3425 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3426 stride_face_dof],
3427 comp,
3428 this->solution_renumbered[i]);
3429
3430 face_dof_values_ptr = this->solution_renumbered.data();
3431 }
3432
3433 constexpr int stride_face_dof_actual =
3434 n_components == 1 ? stride_face_dof : 1;
3435
3436 // loop over quadrature batches qb
3437 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3438
3439 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3440 {
3443
3444 if (evaluation_flags & EvaluationFlags::gradients)
3445 {
3446 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3447 is_linear ?
3449 dim - 1,
3452 2,
3453 stride_face_dof_actual>(face_dof_values_ptr,
3454 this->unit_point_faces_ptr[qb]) :
3456 dim - 1,
3459 2,
3460 false,
3461 stride_face_dof_actual>(this->shapes_faces.data() +
3462 qb * n_shapes,
3463 n_shapes,
3464 face_dof_values_ptr);
3465
3466 value = interpolated_value[dim - 1];
3467 // reorder derivative from tangential/normal derivatives into tensor
3468 // in physical coordinates
3469 if (this->current_face_number / 2 == 0)
3470 {
3471 gradient[0] = interpolated_value[dim];
3472 if (dim > 1)
3473 gradient[1] = interpolated_value[0];
3474 if (dim > 2)
3475 gradient[2] = interpolated_value[1];
3476 }
3477 else if (this->current_face_number / 2 == 1)
3478 {
3479 if (dim > 1)
3480 gradient[1] = interpolated_value[dim];
3481 if (dim == 3)
3482 {
3483 gradient[0] = interpolated_value[1];
3484 gradient[2] = interpolated_value[0];
3485 }
3486 else if (dim == 2)
3487 gradient[0] = interpolated_value[0];
3488 else
3490 }
3491 else if (this->current_face_number / 2 == 2)
3492 {
3493 if (dim > 2)
3494 {
3495 gradient[0] = interpolated_value[0];
3496 gradient[1] = interpolated_value[1];
3497 gradient[2] = interpolated_value[dim];
3498 }
3499 else
3501 }
3502 else
3504 }
3505 else
3506 {
3507 value = is_linear ?
3509 dim - 1,
3512 stride_face_dof_actual>(face_dof_values_ptr,
3513 this->unit_point_faces_ptr[qb]) :
3515 dim - 1,
3518 false,
3519 stride_face_dof_actual>(this->shapes_faces.data() +
3520 qb * n_shapes,
3521 n_shapes,
3522 face_dof_values_ptr);
3523 }
3524
3525 if (evaluation_flags & EvaluationFlags::values)
3526 {
3527 for (unsigned int v = 0, offset = qb * stride;
3528 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3529 ++v, ++offset)
3530 ETT::set_value(value, v, this->values[offset]);
3531 }
3532 if (evaluation_flags & EvaluationFlags::gradients)
3533 {
3534 Assert(this->update_flags & update_gradients ||
3535 this->update_flags & update_inverse_jacobians,
3537
3538 for (unsigned int v = 0, offset = qb * stride;
3539 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3540 ++v, ++offset)
3541 {
3542 unit_gradient_type unit_gradient;
3543 ETT::set_gradient(gradient, v, unit_gradient);
3544 this->gradients[offset] =
3545 this->cell_type <=
3548 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3550 this
3551 ->inverse_jacobian_ptr[this->cell_type <=
3553 GeometryType::affine ?
3554 0 :
3555 offset]
3556 .transpose(),
3557 unit_gradient);
3558 }
3559 }
3560 }
3561}
3562
3563
3564
3565template <int n_components_, int dim, int spacedim, typename Number>
3566template <int stride_face_dof>
3567void
3569 ScalarNumber *face_dof_values,
3570 const EvaluationFlags::EvaluationFlags &integration_flags,
3571 const bool sum_into_values)
3572{
3573 if (this->use_linear_path)
3574 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3575 integration_flags,
3576 sum_into_values);
3577 else
3578 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3579 integration_flags,
3580 sum_into_values);
3581}
3582
3583
3584
3585template <int n_components_, int dim, int spacedim, typename Number>
3586template <bool do_JxW, bool is_linear, int stride_face_dof>
3587inline void
3590 ScalarNumber *face_dof_values,
3591 const EvaluationFlags::EvaluationFlags &integration_flags,
3592 const bool sum_into_values)
3593{
3594 // zero out lanes of incomplete last quadrature point batch
3595 if constexpr (stride == 1)
3596 if (const unsigned int n_filled_lanes =
3597 this->n_q_points_scalar & (n_lanes_internal - 1);
3598 n_filled_lanes > 0)
3599 {
3600 if (integration_flags & EvaluationFlags::values)
3601 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3602 ETT::set_zero_value(this->values.back(), v);
3603 if (integration_flags & EvaluationFlags::gradients)
3604 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3605 ETT::set_zero_gradient(this->gradients.back(), v);
3606 }
3607
3608 std::array<vectorized_value_type,
3609 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3610 solution_values_vectorized_linear = {};
3611
3612 // loop over quadrature batches qb
3613 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3614
3615 const bool cartesian_cell =
3617 const bool affine_cell =
3619 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3620 {
3621 vectorized_value_type value = {};
3623
3624 if (integration_flags & EvaluationFlags::values)
3625 for (unsigned int v = 0, offset = qb * stride;
3626 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3627 ++v, ++offset)
3628 ETT::get_value(value,
3629 v,
3630 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3631 this->values[offset]);
3632
3633 if (integration_flags & EvaluationFlags::gradients)
3634 for (unsigned int v = 0, offset = qb * stride;
3635 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3636 ++v, ++offset)
3637 {
3638 const auto grad_w =
3639 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3640 this->gradients[offset];
3641 ETT::get_gradient(
3642 gradient,
3643 v,
3644 cartesian_cell ?
3645 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3646 grad_w) :
3648 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3649 grad_w));
3650 }
3651
3652 if (integration_flags & EvaluationFlags::gradients)
3653 {
3654 std::array<vectorized_value_type, 2> value_face = {};
3655 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3656
3657 value_face[0] = value;
3658 // fill derivative in physical coordinates into tangential/normal
3659 // derivatives
3660 if (this->current_face_number / 2 == 0)
3661 {
3662 value_face[1] = gradient[0];
3663 if (dim > 1)
3664 gradient_in_face[0] = gradient[1];
3665 if (dim > 2)
3666 gradient_in_face[1] = gradient[2];
3667 }
3668 else if (this->current_face_number / 2 == 1)
3669 {
3670 if (dim > 1)
3671 value_face[1] = gradient[1];
3672 if (dim == 3)
3673 {
3674 gradient_in_face[0] = gradient[2];
3675 gradient_in_face[1] = gradient[0];
3676 }
3677 else if (dim == 2)
3678 gradient_in_face[0] = gradient[0];
3679 else
3681 }
3682 else if (this->current_face_number / 2 == 2)
3683 {
3684 if (dim > 2)
3685 {
3686 value_face[1] = gradient[2];
3687 gradient_in_face[0] = gradient[0];
3688 gradient_in_face[1] = gradient[1];
3689 }
3690 else
3692 }
3693 else
3695
3697 is_linear,
3698 dim - 1,
3701 2>(this->shapes_faces.data() + qb * n_shapes,
3702 n_shapes,
3703 value_face.data(),
3704 gradient_in_face,
3705 is_linear ? solution_values_vectorized_linear.data() :
3706 this->solution_renumbered_vectorized.data(),
3707 this->unit_point_faces_ptr[qb],
3708 qb != 0);
3709 }
3710 else
3712 dim - 1,
3715 this->shapes_faces.data() + qb * n_shapes,
3716 n_shapes,
3717 value,
3718 is_linear ? solution_values_vectorized_linear.data() :
3719 this->solution_renumbered_vectorized.data(),
3720 this->unit_point_faces_ptr[qb],
3721 qb != 0);
3722 }
3723
3724 const unsigned int dofs_per_comp_face =
3725 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3726
3727 for (unsigned int comp = 0; comp < n_components; ++comp)
3728 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3729 if (sum_into_values)
3730 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3731 stride_face_dof] +=
3732 ETT::sum_value(comp,
3733 is_linear ?
3734 *(solution_values_vectorized_linear.data() + i) :
3735 this->solution_renumbered_vectorized[i]);
3736 else
3737 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3738 ETT::sum_value(comp,
3739 is_linear ?
3740 *(solution_values_vectorized_linear.data() + i) :
3741 this->solution_renumbered_vectorized[i]);
3742}
3743
3744
3745
3746template <int n_components_, int dim, int spacedim, typename Number>
3749 const unsigned int point_index) const
3750{
3751 AssertIndexRange(point_index, this->n_q_points);
3752 Assert(this->normal_ptr != nullptr,
3754 ExcFEPointEvaluationAccessToUninitializedMappingField(
3755 "update_normal_vectors"));
3756 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3757 {
3759 for (unsigned int d = 0; d < dim; ++d)
3760 normal[d] =
3761 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3762 0);
3763 if (this->is_interior)
3764 return normal;
3765 else
3766 return -normal;
3767 }
3768 else
3769 {
3770 if (this->is_interior)
3771 return this->normal_ptr[point_index];
3772 else
3773 return -(this->normal_ptr[point_index]);
3774 }
3775}
3776
3778
3779#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
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)
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)