deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
evaluation_kernels_hanging_nodes.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2021 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_evaluation_kernels_hanging_nodes_h
15#define dealii_matrix_free_evaluation_kernels_hanging_nodes_h
16
17#include <deal.II/base/config.h>
18
21
24
25
27
28#ifdef DEBUG
29# define DEAL_II_ALWAYS_INLINE_RELEASE
30#else
31# define DEAL_II_ALWAYS_INLINE_RELEASE DEAL_II_ALWAYS_INLINE
32#endif
33
34
35
36namespace internal
37{
48
49
50
56 enum class HelperType
57 {
68 };
69
70
71
77 {
81 index,
85 group,
90 mask,
94 sorted
95 };
96
97
98
104 int dim,
105 int fe_degree,
106 typename Number>
108
109
110
116 template <int dim, int fe_degree, typename Number>
119 dim,
120 fe_degree,
121 Number>
122 {
123 private:
124 template <int structdim,
125 unsigned int direction,
126 bool transpose,
127 typename Number2>
128 static void
129 interpolate(const unsigned int offset,
130 const unsigned int outer_stride,
131 const unsigned int given_degree,
132 const Number mask_weight,
133 const Number mask_write,
134 const Number2 *DEAL_II_RESTRICT weights,
135 Number *DEAL_II_RESTRICT values)
136 {
137 static constexpr unsigned int max_n_points_1D = 40;
138
139 static_assert(structdim == 1 || structdim == 2,
140 "Only 1D and 2d interpolation implemented");
141 Number temp[fe_degree != -1 ? fe_degree + 1 : max_n_points_1D];
142
143 const unsigned int points =
144 (fe_degree != -1 ? fe_degree : given_degree) + 1;
145
146 AssertIndexRange(points, max_n_points_1D);
147
148 const unsigned int stride = Utilities::pow(points, direction);
149
150 const unsigned int end_of_outer_loop = structdim == 1 ? 2 : points - 1;
151 for (unsigned int g = 1; g < end_of_outer_loop; ++g)
152 {
153 const unsigned int my_offset =
154 offset + (structdim > 1 ? g * outer_stride : 0);
155
156 // extract values to interpolate, setting the first variable to zero
157 // to avoid compile warnings about possibly uninitialized variables
158 temp[0] = 0;
159 for (unsigned int k = 0; k < points; ++k)
160 temp[k] = values[my_offset + k * stride];
161
162 // perform interpolation point by point and write back
163 for (unsigned int k = 0; k < points / 2; ++k)
164 {
165 const unsigned int kmirror = points - 1 - k;
166 Number sum0 = Number(), sum1 = Number(), sum2 = Number(),
167 sum3 = Number();
168 for (unsigned int h = 0; h < points; ++h)
169 {
170 const unsigned int hmirror = points - 1 - h;
171 // load from both sides of the interpolation matrix to
172 // reflect symmetry between the two subfaces along that
173 // direction
174 const Number w0 = weights[(transpose ? 1 : points) * kmirror +
175 (transpose ? points : 1) * hmirror];
176 const Number w1 = weights[(transpose ? 1 : points) * k +
177 (transpose ? points : 1) * h];
178 sum0 += temp[h] * w0;
179 sum1 += temp[h] * w1;
180 sum2 += temp[hmirror] * w1;
181 sum3 += temp[hmirror] * w0;
182 }
183 values[my_offset + k * stride] =
184 temp[k] +
185 mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
186 values[my_offset + kmirror * stride] =
187 temp[kmirror] +
188 mask_write *
189 (sum2 + mask_weight * (sum3 - sum2) - temp[kmirror]);
190 }
191
192 // cleanup case
193 if (points % 2)
194 {
195 const unsigned int k = points / 2;
196 Number sum0 = temp[k] * weights[(transpose ? 1 : points) * k +
197 (transpose ? points : 1) * k],
198 sum1 = sum0;
199 for (unsigned int h = 0; h < points / 2; ++h)
200 {
201 const unsigned int hmirror = points - 1 - h;
202 const Number w0 = weights[(transpose ? 1 : points) * k +
203 (transpose ? points : 1) * hmirror];
204 const Number w1 = weights[(transpose ? 1 : points) * k +
205 (transpose ? points : 1) * h];
206 sum0 += temp[h] * w0;
207 sum0 += temp[hmirror] * w1;
208 sum1 += temp[h] * w1;
209 sum1 += temp[hmirror] * w0;
210 }
211 values[my_offset + k * stride] =
212 temp[k] +
213 mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
214 }
215 }
216 }
217
218 public:
219 template <bool transpose, typename Number2>
220 static void
222 const unsigned int n_components,
225 Number::size()> &constraint_mask,
226 Number *values)
227 {
228 const unsigned int given_degree =
229 fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
230
231 const Number2 *DEAL_II_RESTRICT weights =
232 shape_info.data.front().subface_interpolation_matrices[0].data();
233
234 const unsigned int points = given_degree + 1;
235 const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
236
237 if (dim == 2)
238 {
239 ::ndarray<Number, 2> mask_weights = {};
240 ::ndarray<Number, 2, 2> mask_write = {};
241 ::ndarray<bool, 2, 2> do_face = {};
242
243 for (unsigned int v = 0; v < Number::size(); ++v)
244 {
245 const auto kind = constraint_mask[v];
246 const bool subcell_x = (kind >> 0) & 1;
247 const bool subcell_y = (kind >> 1) & 1;
248 const bool face_x = (kind >> 3) & 1;
249 const bool face_y = (kind >> 4) & 1;
250
251 if (face_y)
252 {
253 const unsigned int side = !subcell_y;
254 mask_write[0][side][v] = 1;
255 do_face[0][side] = true;
256 mask_weights[0][v] = subcell_x;
257 }
258
259 if (face_x)
260 {
261 const unsigned int side = !subcell_x;
262 mask_write[1][side][v] = 1;
263 do_face[1][side] = true;
264 mask_weights[1][v] = subcell_y;
265 }
266 }
267
268 // x direction
269 {
270 const std::array<unsigned int, 2> offsets = {
271 {0, (points - 1) * points}};
272 for (unsigned int c = 0; c < n_components; ++c)
273 for (unsigned int face = 0; face < 2; ++face)
274 if (do_face[0][face])
275 interpolate<1, 0, transpose>(offsets[face],
276 0,
277 given_degree,
278 mask_weights[0],
279 mask_write[0][face],
280 weights,
281 values + c * n_dofs);
282 }
283
284 // y direction
285 {
286 const std::array<unsigned int, 2> offsets = {{0, points - 1}};
287 for (unsigned int c = 0; c < n_components; ++c)
288 for (unsigned int face = 0; face < 2; ++face)
289 if (do_face[1][face])
290 interpolate<1, 1, transpose>(offsets[face],
291 0,
292 given_degree,
293 mask_weights[1],
294 mask_write[1][face],
295 weights,
296 values + c * n_dofs);
297 }
298 }
299 else if (dim == 3)
300 {
301 const unsigned int p0 = 0;
302 const unsigned int p1 = points - 1;
303 const unsigned int p2 = points * points - points;
304 const unsigned int p3 = points * points - 1;
305 const unsigned int p4 = points * points * points - points * points;
306 const unsigned int p5 =
307 points * points * points - points * points + points - 1;
308 const unsigned int p6 = points * points * points - points;
309
310 ::ndarray<bool, 3, 4> process_edge = {};
311 ::ndarray<bool, 3, 4> process_face = {};
312 ::ndarray<Number, 3, 4> mask_edge = {};
313 ::ndarray<Number, 3, 4> mask_face = {};
314 ::ndarray<Number, 3> mask_weights = {};
315
316 for (unsigned int v = 0; v < Number::size(); ++v)
317 {
318 const auto kind = constraint_mask[v];
319
320 const bool subcell_x = (kind >> 0) & 1;
321 const bool subcell_y = (kind >> 1) & 1;
322 const bool subcell_z = (kind >> 2) & 1;
323 const bool face_x = ((kind >> 3) & 1) ? (kind >> 5) & 1 : 0;
324 const bool face_y = ((kind >> 3) & 1) ? (kind >> 6) & 1 : 0;
325 const bool face_z = ((kind >> 3) & 1) ? (kind >> 7) & 1 : 0;
326 const bool edge_x = ((kind >> 4) & 1) ? (kind >> 5) & 1 : 0;
327 const bool edge_y = ((kind >> 4) & 1) ? (kind >> 6) & 1 : 0;
328 const bool edge_z = ((kind >> 4) & 1) ? (kind >> 7) & 1 : 0;
329
330 if (subcell_x)
331 mask_weights[0][v] = 1;
332 if (subcell_y)
333 mask_weights[1][v] = 1;
334 if (subcell_z)
335 mask_weights[2][v] = 1;
336
337 if (face_x)
338 {
339 const unsigned int side = !subcell_x;
340
341 mask_face[1][side][v] = process_face[1][side] = true;
342 mask_edge[1][side][v] = process_edge[1][side] = true;
343 mask_edge[1][2 + side][v] = process_edge[1][2 + side] = true;
344 mask_face[2][side][v] = process_face[2][side] = true;
345 mask_edge[2][side][v] = process_edge[2][side] = true;
346 mask_edge[2][2 + side][v] = process_edge[2][2 + side] = true;
347 }
348 if (face_y)
349 {
350 const unsigned int side = !subcell_y;
351
352 mask_face[0][side][v] = process_face[0][side] = true;
353 mask_edge[0][side][v] = process_edge[0][side] = true;
354 mask_edge[0][2 + side][v] = process_edge[0][2 + side] = true;
355 mask_face[2][2 + side][v] = process_face[2][2 + side] = true;
356 mask_edge[2][2 * side][v] = process_edge[2][2 * side] = true;
357 mask_edge[2][2 * side + 1][v] =
358 process_edge[2][2 * side + 1] = true;
359 }
360 if (face_z)
361 {
362 const unsigned int side = !subcell_z;
363
364 mask_face[0][2 + side][v] = process_face[0][2 + side] = true;
365 mask_edge[0][2 * side][v] = process_edge[0][2 * side] = true;
366 mask_edge[0][2 * side + 1][v] =
367 process_edge[0][2 * side + 1] = true;
368 mask_face[1][2 + side][v] = process_face[1][2 + side] = true;
369 mask_edge[1][2 * side][v] = process_edge[1][2 * side] = true;
370 mask_edge[1][2 * side + 1][v] =
371 process_edge[1][2 * side + 1] = true;
372 }
373 if (edge_x)
374 {
375 const unsigned int index = (!subcell_z) * 2 + (!subcell_y);
376 mask_edge[0][index][v] = process_edge[0][index] = true;
377 }
378 if (edge_y)
379 {
380 const unsigned int index = (!subcell_z) * 2 + (!subcell_x);
381 mask_edge[1][index][v] = process_edge[1][index] = true;
382 }
383 if (edge_z)
384 {
385 const unsigned int index = (!subcell_y) * 2 + (!subcell_x);
386 mask_edge[2][index][v] = process_edge[2][index] = true;
387 }
388 }
389
390 // direction 0:
391 if (given_degree > 1)
392 {
393 const std::array<unsigned int, 4> face_offsets = {
394 {p0, p2, p0, p4}};
395 const std::array<unsigned int, 2> outer_strides = {
396 {points * points, points}};
397 for (unsigned int c = 0; c < n_components; ++c)
398 for (unsigned int face = 0; face < 4; ++face)
399 if (process_face[0][face])
400 interpolate<2, 0, transpose>(face_offsets[face],
401 outer_strides[face / 2],
402 given_degree,
403 mask_weights[0],
404 mask_face[0][face],
405 weights,
406 values + c * n_dofs);
407 }
408 {
409 const std::array<unsigned int, 4> edge_offsets = {{p0, p2, p4, p6}};
410 for (unsigned int c = 0; c < n_components; ++c)
411 for (unsigned int edge = 0; edge < 4; ++edge)
412 if (process_edge[0][edge])
413 interpolate<1, 0, transpose>(edge_offsets[edge],
414 0,
415 given_degree,
416 mask_weights[0],
417 mask_edge[0][edge],
418 weights,
419 values + c * n_dofs);
420 }
421
422 // direction 1:
423 if (given_degree > 1)
424 {
425 const std::array<unsigned int, 4> face_offsets = {
426 {p0, p1, p0, p4}};
427 const std::array<unsigned int, 2> outer_strides = {
428 {points * points, 1}};
429 for (unsigned int c = 0; c < n_components; ++c)
430 for (unsigned int face = 0; face < 4; ++face)
431 if (process_face[1][face])
432 interpolate<2, 1, transpose>(face_offsets[face],
433 outer_strides[face / 2],
434 given_degree,
435 mask_weights[1],
436 mask_face[1][face],
437 weights,
438 values + c * n_dofs);
439 }
440
441 {
442 const std::array<unsigned int, 4> edge_offsets = {{p0, p1, p4, p5}};
443 for (unsigned int c = 0; c < n_components; ++c)
444 for (unsigned int edge = 0; edge < 4; ++edge)
445 if (process_edge[1][edge])
446 interpolate<1, 1, transpose>(edge_offsets[edge],
447 0,
448 given_degree,
449 mask_weights[1],
450 mask_edge[1][edge],
451 weights,
452 values + c * n_dofs);
453 }
454
455 // direction 2:
456 if (given_degree > 1)
457 {
458 const std::array<unsigned int, 4> face_offsets = {
459 {p0, p1, p0, p2}};
460 const std::array<unsigned int, 2> outer_strides = {{points, 1}};
461 for (unsigned int c = 0; c < n_components; ++c)
462 for (unsigned int face = 0; face < 4; ++face)
463 if (process_face[2][face])
464 interpolate<2, 2, transpose>(face_offsets[face],
465 outer_strides[face / 2],
466 given_degree,
467 mask_weights[2],
468 mask_face[2][face],
469 weights,
470 values + c * n_dofs);
471 }
472
473 {
474 const std::array<unsigned int, 4> edge_offsets = {{p0, p1, p2, p3}};
475 for (unsigned int c = 0; c < n_components; ++c)
476 for (unsigned int edge = 0; edge < 4; ++edge)
477 if (process_edge[2][edge])
478 interpolate<1, 2, transpose>(edge_offsets[edge],
479 0,
480 given_degree,
481 mask_weights[2],
482 mask_edge[2][edge],
483 weights,
484 values + c * n_dofs);
485 }
486 }
487 else
488 {
490 }
491 }
492 };
493
494 template <typename T1, VectorizationTypes VT>
495 struct Trait;
496
497 template <typename T1>
499 {
500 using value_type = typename T1::value_type;
501 using index_type = unsigned int;
503
504 template <typename T>
505 static inline const std::array<AlignedVector<interpolation_type>, 2> &
506 get_interpolation_matrix(const T &shape_info)
507 {
508 return shape_info.data.front().subface_interpolation_matrices_scalar;
509 }
510
511 static inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
513 T1::size()> mask,
515 T1::size()> mask_new,
516 const unsigned int v)
517 {
518 (void)mask;
519 (void)mask_new;
520 return v;
521 }
522
523 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
524 do_break(unsigned int v,
526 {
527 (void)v;
528 (void)kind;
529 return false;
530 }
531
532 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
533 do_continue(unsigned int v,
535 {
536 (void)v;
537 return kind ==
539 }
540
545 T1::size()> mask)
546 {
547 return mask;
548 }
549
550 static inline DEAL_II_ALWAYS_INLINE_RELEASE typename T1::value_type
551 get_value(const typename T1::value_type &value, const index_type &i)
552 {
553 (void)i;
554 return value;
555 }
556
557 static inline DEAL_II_ALWAYS_INLINE_RELEASE typename T1::value_type
558 get_value(const T1 &value, const index_type &i)
559 {
560 return value[i];
561 }
562
563 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
564 set_value(T1 &result,
565 const typename T1::value_type &value,
566 const index_type &i)
567 {
568 result[i] = value;
569 }
570 };
571
572 template <typename T1>
574 {
575 using value_type = T1;
576 using index_type = std::pair<T1, T1>;
578
579 template <typename T>
580 static inline const std::array<AlignedVector<T1>, 2> &
581 get_interpolation_matrix(const T &shape_info)
582 {
583 return shape_info.data.front().subface_interpolation_matrices;
584 }
585
586 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
587 do_break(unsigned int v,
589 {
590 (void)v;
591 (void)kind;
592 return false;
593 }
594
595 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
596 do_continue(unsigned int v,
598 {
599 (void)v;
600 return kind ==
602 }
603
604 static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type
606 T1::size()> mask,
608 T1::size()> mask_new,
609 const unsigned int v)
610 {
611 (void)mask;
612 (void)mask_new;
613 T1 result = 0.0;
614 result[v] = 1.0;
615 return {result, T1(1.0) - result};
616 }
617
622 T1::size()> mask)
623 {
624 return mask;
625 }
626
627 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
628 get_value(const T1 &value, const index_type &)
629 {
630 return value;
631 }
632
633 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
634 set_value(T1 &result, const T1 &value, const index_type &i)
635 {
636 result = result * i.second + value * i.first;
637 }
638 };
639
640 template <typename T1>
642 {
643 using value_type = T1;
644 using index_type = std::pair<T1, T1>;
646
647 template <typename T>
648 static inline const std::array<AlignedVector<T1>, 2> &
649 get_interpolation_matrix(const T &shape_info)
650 {
651 return shape_info.data.front().subface_interpolation_matrices;
652 }
653
654 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
655 do_break(unsigned int v,
657 {
658 (void)v;
659 return kind ==
661 }
662
663 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
664 do_continue(unsigned int v,
666 {
667 (void)v;
668 return kind ==
670 }
671
672 static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type
674 T1::size()> mask,
676 T1::size()> mask_new,
677 const unsigned int v)
678 {
679 T1 result;
680
681 for (unsigned int i = 0; i < T1::size(); ++i)
682 result[i] = mask_new[v] == mask[i];
683
684 return {result, T1(1.0) - result};
685 }
686
691 T1::size()> mask)
692 {
693 auto new_mask = mask;
694
695 std::sort(new_mask.begin(), new_mask.end());
696 std::fill(std::unique(new_mask.begin(), new_mask.end()),
697 new_mask.end(),
699
700 return new_mask;
701 }
702
703 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
704 get_value(const T1 &value, const index_type &)
705 {
706 return value;
707 }
708
709 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
710 set_value(T1 &result, const T1 &value, const index_type &i)
711 {
712 result = result * i.second + value * i.first;
713 }
714 };
715
716 template <typename T1>
718 {
719 using value_type = T1;
720 using index_type = T1;
722
723 template <typename T>
724 static inline const std::array<AlignedVector<T1>, 2> &
725 get_interpolation_matrix(const T &shape_info)
726 {
727 return shape_info.data.front().subface_interpolation_matrices;
728 }
729
730 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
731 do_break(unsigned int v,
733 {
734 (void)kind;
735 return v > 0;
736 }
737
738 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
739 do_continue(unsigned int v,
741 {
742 (void)kind;
743
745
746 return v > 0; // should not be called
747 }
748
749 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
751 T1::size()> mask,
753 T1::size()> mask_new,
754 const unsigned int v)
755 {
756 (void)mask;
757 (void)mask_new;
758 (void)v;
759 return 1.0; // return something since not used
760 }
761
766 T1::size()> mask)
767 {
768 return mask;
769 }
770
771 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
772 get_value(const T1 &value, const index_type &)
773 {
774 return value;
775 }
776
777 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
778 set_value(T1 &result, const T1 &value, const index_type &i)
779 {
780 (void)i;
781 result = value;
782 }
783 };
784
785
786
797 template <typename T,
798 typename Number,
799 VectorizationTypes VectorizationType,
800 int fe_degree,
801 bool transpose>
803 {
804 public:
807 const T &t,
808 const unsigned int &given_degree,
809 const bool &type_x,
810 const bool &type_y,
811 const bool &type_z,
813 const std::array<
817 Number *values)
818 : t(t)
820 , type_x(type_x)
821 , type_y(type_y)
822 , type_z(type_z)
823 , v(v)
825 , values(values)
826 {}
827
828 template <unsigned int direction, unsigned int d, bool skip_borders>
829 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
831 const unsigned int dof_offset,
832 const unsigned int given_degree,
835 *DEAL_II_RESTRICT weight,
836 Number *DEAL_II_RESTRICT values)
837 {
838 static constexpr unsigned int max_n_points_1D = 40;
839
841 temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
842
843 const unsigned int points =
844 (fe_degree != -1 ? fe_degree : given_degree) + 1;
845
846 AssertIndexRange(given_degree, max_n_points_1D);
847
848 const unsigned int stride = fe_degree != -1 ?
849 Utilities::pow(fe_degree + 1, direction) :
850 Utilities::pow(given_degree + 1, direction);
851
852 // direction side0 side1 side2
853 // 0 - p^2 p
854 // 1 p^2 - 1
855 // 2 p - 1
856 const unsigned int stride2 =
857 ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ?
858 (points * points) :
859 (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points :
860 1);
861
862 for (unsigned int g = (skip_borders ? 1 : 0);
863 g < points - (skip_borders ? 1 : 0);
864 ++g)
865 {
866 // copy result back
867 for (unsigned int k = 0; k < points; ++k)
869 values[dof_offset + k * stride + stride2 * g], v);
870
871 // perform interpolation point by point
872 for (unsigned int k = 0; k < points; ++k)
873 {
875 weight[(transpose ? 1 : points) * k], v) *
876 temp[0];
877 for (unsigned int h = 1; h < points; ++h)
879 weight[(transpose ? 1 : points) * k +
880 (transpose ? points : 1) * h],
881 v) *
882 temp[h];
884 values[dof_offset + k * stride + stride2 * g], sum, v);
885 }
886 }
887 }
888
889 template <unsigned int direction>
890 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
892 const unsigned int p,
893 const unsigned int given_degree,
896 *DEAL_II_RESTRICT weight,
897 Number *DEAL_II_RESTRICT values)
898 {
899 static constexpr unsigned int max_n_points_1D = 40;
900
902 temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
903
904 const unsigned int points =
905 (fe_degree != -1 ? fe_degree : given_degree) + 1;
906
907 AssertIndexRange(given_degree, max_n_points_1D);
908
909 const unsigned int stride = fe_degree != -1 ?
910 Utilities::pow(fe_degree + 1, direction) :
911 Utilities::pow(given_degree + 1, direction);
912
913 // copy result back
914 for (unsigned int k = 0; k < points; ++k)
915 temp[k] =
917 v);
918
919 // perform interpolation point by point
920 for (unsigned int k = 0; k < points; ++k)
921 {
923 weight[(transpose ? 1 : points) * k], v) *
924 temp[0];
925 for (unsigned int h = 1; h < points; ++h)
927 weight[(transpose ? 1 : points) * k +
928 (transpose ? points : 1) * h],
929 v) *
930 temp[h];
932 sum,
933 v);
934 }
935 }
936
937 template <bool do_x, bool do_y, bool do_z>
940 {
941 if (do_x)
942 interpolate_3D_edge<0>(t.line(0, type_y, type_z),
944 v,
946 values);
947
948 if (do_y)
949 interpolate_3D_edge<1>(t.line(1, type_x, type_z),
951 v,
953 values);
954
955 if (do_z)
956 interpolate_3D_edge<2>(t.line(2, type_x, type_y),
958 v,
960 values);
961 }
962
963 template <bool do_x, bool do_y, bool do_z>
966 {
967 static_assert((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) ||
968 (!do_x && !do_y && do_z),
969 "Only one face can be chosen.");
970
971 static const unsigned int direction = do_x ? 0 : (do_y ? 1 : 2);
972 const bool type = do_x ? type_x : (do_y ? type_y : type_z);
973
974 if (!do_x)
975 interpolate_3D_face<0, direction, false>(
976 t.face(direction, type),
978 v,
980 values);
981
982 if (!do_y)
983 interpolate_3D_face<1, direction, false>(
984 t.face(direction, type),
986 v,
988 values);
989
990 if (!do_z)
991 interpolate_3D_face<2, direction, false>(
992 t.face(direction, type),
994 v,
996 values);
997 }
998
999 template <bool do_x, bool do_y, bool do_z>
1002 {
1003 static_assert(((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) ||
1004 (!do_x && !do_y && do_z)) == false,
1005 "Only one face can be chosen.");
1006
1007 // direction 0
1008 {
1009 const auto inpterolation_matrix =
1011
1012 // faces
1013 if (do_y && given_degree > 1)
1014 interpolate_3D_face<0, 1, true>(
1015 t.face(1, type_y), given_degree, v, inpterolation_matrix, values);
1016
1017 if (do_z && given_degree > 1)
1018 interpolate_3D_face<0, 2, true>(
1019 t.face(2, type_z), given_degree, v, inpterolation_matrix, values);
1020
1021 // direction 0 -> edges
1022 interpolate_3D_edge<0>((do_x && do_y && !do_z) ?
1023 (t.lines_plane(0, type_x, type_y, 0)) :
1024 ((do_x && !do_y && do_z) ?
1025 (t.lines_plane(1, type_x, type_z, 0)) :
1026 (t.lines(0, type_y, type_z, 0))),
1028 v,
1029 inpterolation_matrix,
1030 values);
1031
1032
1033 interpolate_3D_edge<0>((do_x && do_y && !do_z) ?
1034 (t.lines_plane(0, type_x, type_y, 1)) :
1035 ((do_x && !do_y && do_z) ?
1036 (t.lines_plane(1, type_x, type_z, 1)) :
1037 (t.lines(0, type_y, type_z, 1))),
1039 v,
1040 inpterolation_matrix,
1041 values);
1042
1043 if (do_y && do_z)
1044 interpolate_3D_edge<0>(t.lines(0, type_y, type_z, 2),
1046 v,
1047 inpterolation_matrix,
1048 values);
1049 }
1050
1051 // direction 1
1052 {
1053 const auto inpterolation_matrix =
1055
1056 // faces
1057 if (do_x && given_degree > 1)
1058 interpolate_3D_face<1, 0, true>(
1059 t.face(0, type_x), given_degree, v, inpterolation_matrix, values);
1060
1061 if (do_z && given_degree > 1)
1062 interpolate_3D_face<1, 2, true>(
1063 t.face(2, type_z), given_degree, v, inpterolation_matrix, values);
1064
1065 // lines
1066 interpolate_3D_edge<1>((do_x && do_y && !do_z) ?
1067 (t.lines_plane(0, type_x, type_y, 2)) :
1068 ((!do_x && do_y && do_z) ?
1069 (t.lines_plane(2, type_y, type_z, 0)) :
1070 (t.lines(1, type_x, type_z, 0))),
1072 v,
1073 inpterolation_matrix,
1074 values);
1075
1076 interpolate_3D_edge<1>((do_x && do_y && !do_z) ?
1077 (t.lines_plane(0, type_x, type_y, 3)) :
1078 ((!do_x && do_y && do_z) ?
1079 (t.lines_plane(2, type_y, type_z, 1)) :
1080 (t.lines(1, type_x, type_z, 1))),
1082 v,
1083 inpterolation_matrix,
1084 values);
1085
1086 if (do_x && do_z)
1087 interpolate_3D_edge<1>(t.lines(1, type_x, type_z, 2),
1089 v,
1090 inpterolation_matrix,
1091 values);
1092 }
1093
1094 // direction 2 -> faces
1095 {
1096 const auto inpterolation_matrix =
1098
1099 if (do_x && given_degree > 1)
1100 interpolate_3D_face<2, 0, true>(
1101 t.face(0, type_x), given_degree, v, inpterolation_matrix, values);
1102
1103 if (do_y && given_degree > 1)
1104 interpolate_3D_face<2, 1, true>(
1105 t.face(1, type_y), given_degree, v, inpterolation_matrix, values);
1106
1107 // direction 2 -> edges
1108 interpolate_3D_edge<2>((do_x && !do_y && do_z) ?
1109 (t.lines_plane(1, type_x, type_z, 2)) :
1110 ((!do_x && do_y && do_z) ?
1111 (t.lines_plane(2, type_y, type_z, 2)) :
1112 (t.lines(2, type_x, type_y, 0))),
1114 v,
1115 inpterolation_matrix,
1116 values);
1117
1118 interpolate_3D_edge<2>((do_x && !do_y && do_z) ?
1119 (t.lines_plane(1, type_x, type_z, 3)) :
1120 ((!do_x && do_y && do_z) ?
1121 (t.lines_plane(2, type_y, type_z, 3)) :
1122 (t.lines(2, type_x, type_y, 1))),
1124 v,
1125 inpterolation_matrix,
1126 values);
1127
1128 if (do_x && do_y)
1129 interpolate_3D_edge<2>(t.lines(2, type_x, type_y, 2),
1131 v,
1132 inpterolation_matrix,
1133 values);
1134 }
1135 }
1136
1137 private:
1138 const T &t;
1139 const unsigned int &given_degree;
1140 const bool &type_x;
1141 const bool &type_y;
1142 const bool &type_z;
1144 const std::array<
1148 Number *values;
1149 };
1150
1155 template <HelperType helper_type,
1156 typename Number,
1157 VectorizationTypes VectorizationType,
1158 int fe_degree,
1159 bool transpose>
1161
1171 template <typename Number,
1172 VectorizationTypes VectorizationType,
1173 int fe_degree,
1174 bool transpose>
1177 Number,
1178 VectorizationType,
1179 fe_degree,
1180 transpose>
1182 FEEvaluationImplHangingNodesScalarEntityInterpolationImpl<
1183 HelperType::dynamic,
1184 Number,
1185 VectorizationType,
1186 fe_degree,
1187 transpose>,
1188 Number,
1189 VectorizationType,
1190 fe_degree,
1191 transpose>
1192 {
1193 public:
1194 // Compiling with gcc 16.1.1 results in an annoying warning for the
1195 // following constructor: <unknown>’ may be used uninitialized
1196 // This may be due to (*this) being passed to the base class constructor.
1197 // We disable this warning.
1201 const unsigned int &given_degree,
1202 const bool &type_x,
1203 const bool &type_y,
1204 const bool &type_z,
1206 const std::array<
1209 2> &interpolation_matrices,
1210 Number *values)
1214 Number,
1215 VectorizationType,
1216 fe_degree,
1217 transpose>,
1218 Number,
1219 VectorizationType,
1220 fe_degree,
1221 transpose>(*this,
1222 given_degree,
1223 type_x,
1224 type_y,
1225 type_z,
1226 v,
1227 interpolation_matrices,
1228 values)
1229 , points(given_degree + 1)
1230 {
1231 static_assert(fe_degree == -1, "Only working for fe_degree = -1.");
1232 }
1234
1235 const unsigned int points;
1236
1237 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1238 line(unsigned int i, unsigned int j, unsigned int k) const
1239 {
1240 return line_array[i][j][k];
1241 }
1242
1243 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1244 face(unsigned int i, unsigned int j) const
1245 {
1246 return face_array[i][j];
1247 }
1248
1249 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1250 lines_plane(unsigned int i,
1251 unsigned int j,
1252 unsigned int k,
1253 unsigned int l) const
1254 {
1255 return lines_plane_array[i][j][k][l];
1256 }
1257
1258 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1259 lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
1260 {
1261 return lines_array[i][j][k][l];
1262 }
1263
1264 private:
1265 const ::ndarray<unsigned int, 3, 2, 2> line_array = {
1266 {{{{{points * points * points - points, points *points - points}},
1267 {{points * points * points - points * points, 0}}}},
1268 {{{{points * points * points - points * points + points - 1,
1269 points - 1}},
1270 {{points * points * points - points * points, 0}}}},
1271 {{{{points * points - 1, points - 1}},
1272 {{points * points - points, 0}}}}}};
1273
1274 const ::ndarray<unsigned int, 3, 2> face_array = {
1275 {{{points - 1, 0}},
1276 {{points * points - points, 0}},
1277 {{points * points * points - points * points, 0}}}};
1278
1279 const ::ndarray<unsigned int, 3, 2, 2, 4> lines_plane_array = {
1280 {{{{{{{points * points - points,
1281 points *points *points - points,
1282 points - 1,
1283 points *points *points - points *points + points - 1}},
1284 {{0,
1285 points *points *points - points *points,
1286 points - 1,
1287 points *points *points - points *points + points - 1}}}},
1288 {{{{points * points - points,
1289 points *points *points - points,
1290 0,
1291 points *points *points - points *points}},
1292 {{0,
1293 points *points *points - points *points,
1294 0,
1295 points *points *points - points *points}}}}}},
1296 {{{{{{points * points * points - points * points,
1297 points *points *points - points,
1298 points - 1,
1299 points *points - 1}},
1300 {{0, points *points - points, points - 1, points *points - 1}}}},
1301 {{{{points * points * points - points * points,
1302 points *points *points - points,
1303 0,
1304 points *points - points}},
1305 {{0, points *points - points, 0, points *points - points}}}}}},
1306 {{{{{{points * points * points - points * points,
1307 points *points *points - points *points + points - 1,
1308 points *points - points,
1309 points *points - 1}},
1310 {{0, points - 1, points *points - points, points *points - 1}}}},
1311 {{{{points * points * points - points * points,
1312 points *points *points - points *points + points - 1,
1313 0,
1314 points - 1}},
1315 {{0, points - 1, 0, points - 1}}}}}}}};
1316
1317 const ::ndarray<unsigned int, 3, 2, 2, 3> lines_array = {
1318 {{{{{{{points * points - points,
1319 points *points *points - points *points,
1320 points *points *points - points}},
1321 {{0, points *points - points, points *points *points - points}}}},
1322 {{{{0,
1323 points *points *points - points *points,
1324 points *points *points - points}},
1325 {{0,
1326 points *points - points,
1327 points *points *points - points *points}}}}}},
1328 {{{{{{points - 1,
1329 points *points *points - points *points,
1330 points *points *points - points *points + points - 1}},
1331 {{0,
1332 points - 1,
1333 points *points *points - points *points + points - 1}}}},
1334 {{{{0,
1335 points *points *points - points *points,
1336 points *points *points - points *points + points - 1}},
1337 {{0, points - 1, points *points *points - points *points}}}}}},
1338 {{{{{{points - 1, points *points - points, points *points - 1}},
1339 {{0, points - 1, points *points - 1}}}},
1340 {{{{0, points *points - points, points *points - 1}},
1341 {{0, points - 1, points *points - points}}}}}}}};
1342 };
1343
1353 template <typename Number,
1354 VectorizationTypes VectorizationType,
1355 int fe_degree,
1356 bool transpose>
1359 Number,
1360 VectorizationType,
1361 fe_degree,
1362 transpose>
1364 FEEvaluationImplHangingNodesScalarEntityInterpolationImpl<
1365 HelperType::constant,
1366 Number,
1367 VectorizationType,
1368 fe_degree,
1369 transpose>,
1370 Number,
1371 VectorizationType,
1372 fe_degree,
1373 transpose>
1374 {
1375 public:
1376 // Compiling with gcc 16.1.1 results in an annoying warning for the
1377 // following constructor: <unknown>’ may be used uninitialized
1378 // This may be due to (*this) being passed to the base class constructor.
1379 // We disable this warning.
1383 const unsigned int &given_degree,
1384 const bool &type_x,
1385 const bool &type_y,
1386 const bool &type_z,
1388 const std::array<
1391 2> &interpolation_matrices,
1392 Number *values)
1396 Number,
1397 VectorizationType,
1398 fe_degree,
1399 transpose>,
1400 Number,
1401 VectorizationType,
1402 fe_degree,
1403 transpose>(*this,
1404 given_degree,
1405 type_x,
1406 type_y,
1407 type_z,
1408 v,
1409 interpolation_matrices,
1410 values)
1411 {
1412 static_assert(fe_degree != -1, "Only working for fe_degree != -1.");
1413 }
1415
1416
1417 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1418 line(unsigned int i, unsigned int j, unsigned int k) const
1419 {
1420 static constexpr unsigned int points = fe_degree + 1;
1421
1422 static constexpr ::ndarray<unsigned int, 3, 2, 2> line_array = {
1423 {{{{{points * points * points - points, points * points - points}},
1424 {{points * points * points - points * points, 0}}}},
1425 {{{{points * points * points - points * points + points - 1,
1426 points - 1}},
1427 {{points * points * points - points * points, 0}}}},
1428 {{{{points * points - 1, points - 1}},
1429 {{points * points - points, 0}}}}}};
1430
1431 return line_array[i][j][k];
1432 }
1433
1434 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1435 face(unsigned int i, unsigned int j) const
1436 {
1437 static constexpr unsigned int points = fe_degree + 1;
1438
1439 static constexpr ::ndarray<unsigned int, 3, 2> face_array = {
1440 {{{points - 1, 0}},
1441 {{points * points - points, 0}},
1442 {{points * points * points - points * points, 0}}}};
1443
1444 return face_array[i][j];
1445 }
1446
1447 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1448 lines_plane(unsigned int i,
1449 unsigned int j,
1450 unsigned int k,
1451 unsigned int l) const
1452 {
1453 static constexpr unsigned int points = fe_degree + 1;
1454
1455 static constexpr ::ndarray<unsigned int, 3, 2, 2, 4>
1456 lines_plane_array = {
1457 {{{{{{{points * points - points,
1458 points * points * points - points,
1459 points - 1,
1460 points * points * points - points * points + points - 1}},
1461 {{0,
1462 points * points * points - points * points,
1463 points - 1,
1464 points * points * points - points * points + points - 1}}}},
1465 {{{{points * points - points,
1466 points * points * points - points,
1467 0,
1468 points * points * points - points * points}},
1469 {{0,
1470 points * points * points - points * points,
1471 0,
1472 points * points * points - points * points}}}}}},
1473 {{{{{{points * points * points - points * points,
1474 points * points * points - points,
1475 points - 1,
1476 points * points - 1}},
1477 {{0,
1478 points * points - points,
1479 points - 1,
1480 points * points - 1}}}},
1481 {{{{points * points * points - points * points,
1482 points * points * points - points,
1483 0,
1484 points * points - points}},
1485 {{0, points * points - points, 0, points * points - points}}}}}},
1486 {{{{{{points * points * points - points * points,
1487 points * points * points - points * points + points - 1,
1488 points * points - points,
1489 points * points - 1}},
1490 {{0,
1491 points - 1,
1492 points * points - points,
1493 points * points - 1}}}},
1494 {{{{points * points * points - points * points,
1495 points * points * points - points * points + points - 1,
1496 0,
1497 points - 1}},
1498 {{0, points - 1, 0, points - 1}}}}}}}};
1499
1500 return lines_plane_array[i][j][k][l];
1501 }
1502
1503 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1504 lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
1505 {
1506 static constexpr unsigned int points = fe_degree + 1;
1507
1508 static constexpr ::ndarray<unsigned int, 3, 2, 2, 3> lines_array = {
1509 {{{{{{{points * points - points,
1510 points * points * points - points * points,
1511 points * points * points - points}},
1512 {{0,
1513 points * points - points,
1514 points * points * points - points}}}},
1515 {{{{0,
1516 points * points * points - points * points,
1517 points * points * points - points}},
1518 {{0,
1519 points * points - points,
1520 points * points * points - points * points}}}}}},
1521 {{{{{{points - 1,
1522 points * points * points - points * points,
1523 points * points * points - points * points + points - 1}},
1524 {{0,
1525 points - 1,
1526 points * points * points - points * points + points - 1}}}},
1527 {{{{0,
1528 points * points * points - points * points,
1529 points * points * points - points * points + points - 1}},
1530 {{0, points - 1, points * points * points - points * points}}}}}},
1531 {{{{{{points - 1, points * points - points, points * points - 1}},
1532 {{0, points - 1, points * points - 1}}}},
1533 {{{{0, points * points - points, points * points - 1}},
1534 {{0, points - 1, points * points - points}}}}}}}};
1535
1536 return lines_array[i][j][k][l];
1537 }
1538 };
1539
1540
1546 template <int dim, int fe_degree, typename Number>
1549 dim,
1550 fe_degree,
1551 Number>
1552 {
1553 public:
1554 static const VectorizationTypes VectorizationType =
1556
1557 private:
1558 template <unsigned int side, bool transpose>
1559 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
1561 const unsigned int given_degree,
1564 *DEAL_II_RESTRICT weight,
1565 Number *DEAL_II_RESTRICT values)
1566 {
1567 static constexpr unsigned int max_n_points_1D = 40;
1568
1570 temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
1571
1572 const unsigned int points =
1573 (fe_degree != -1 ? fe_degree : given_degree) + 1;
1574
1575 AssertIndexRange(given_degree, max_n_points_1D);
1576
1577 const unsigned int d = side / 2; // direction
1578 const unsigned int s = side % 2; // left or right surface
1579
1580 const unsigned int offset = ::Utilities::pow(points, d + 1);
1581 const unsigned int stride =
1582 (s == 0 ? 0 : (points - 1)) * ::Utilities::pow(points, d);
1583
1584 const unsigned int r1 = ::Utilities::pow(points, dim - d - 1);
1585 const unsigned int r2 = ::Utilities::pow(points, d);
1586
1587 // copy result back
1588 for (unsigned int i = 0, k = 0; i < r1; ++i)
1589 for (unsigned int j = 0; j < r2; ++j, ++k)
1591 values[i * offset + stride + j], v);
1592
1593 // perform interpolation point by point (note: r1 * r2 ==
1594 // points^(dim-1))
1595 for (unsigned int i = 0, k = 0; i < r1; ++i)
1596 for (unsigned int j = 0; j < r2; ++j, ++k)
1597 {
1599 for (unsigned int h = 0; h < points; ++h)
1601 weight[(transpose ? 1 : points) * k +
1602 (transpose ? points : 1) * h],
1603 v) *
1604 temp[h];
1606 values[i * offset + stride + j], sum, v);
1607 }
1608 }
1609
1610 public:
1611 template <bool transpose, typename Number2>
1612 static void
1614 const unsigned int n_desired_components,
1617 Number::size()> &constraint_mask,
1618 Number *values)
1619 {
1620 const unsigned int given_degree =
1621 fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
1622
1623 const auto &interpolation_matrices =
1625
1626 const auto constraint_mask_sorted =
1628
1629 for (unsigned int c = 0; c < n_desired_components; ++c)
1630 {
1631 for (unsigned int v = 0; v < Number::size(); ++v)
1632 {
1633 const auto mask = constraint_mask_sorted[v];
1634
1636 break;
1637
1639 continue;
1640
1641 const auto vv =
1643 constraint_mask_sorted,
1644 v);
1645
1646 if (dim == 2) // 2d: only faces
1647 {
1648 const bool subcell_x = (mask >> 0) & 1;
1649 const bool subcell_y = (mask >> 1) & 1;
1650 const bool face_x = (mask >> 3) & 1;
1651 const bool face_y = (mask >> 4) & 1;
1652
1653 // direction 0:
1654 if (face_y)
1655 {
1656 const auto *weights =
1657 interpolation_matrices[!subcell_x].data();
1658
1659 if (subcell_y)
1660 interpolate_2D<2, transpose>(given_degree,
1661 vv,
1662 weights,
1663 values); // face 2
1664 else
1665 interpolate_2D<3, transpose>(given_degree,
1666 vv,
1667 weights,
1668 values); // face 3
1669 }
1670
1671 // direction 1:
1672 if (face_x)
1673 {
1674 const auto *weights =
1675 interpolation_matrices[!subcell_y].data();
1676
1677 if (subcell_x)
1678 interpolate_2D<0, transpose>(given_degree,
1679 vv,
1680 weights,
1681 values); // face 0
1682 else
1683 interpolate_2D<1, transpose>(given_degree,
1684 vv,
1685 weights,
1686 values); // face 1
1687 }
1688 }
1689 else if (dim == 3) // 3d faces and edges
1690 {
1691 const bool type_x = (mask >> 0) & 1;
1692 const bool type_y = (mask >> 1) & 1;
1693 const bool type_z = (mask >> 2) & 1;
1694
1695 const auto flag_0 = (mask >> 3) & 3;
1696 const auto flag_1 = (mask >> 5) & 7;
1697 const auto faces = (flag_0 & 0b01) ? flag_1 : 0;
1698 const auto edges = (flag_0 & 0b10) ? flag_1 : 0;
1699
1701 fe_degree == -1 ? HelperType::dynamic :
1703 Number,
1704 VectorizationType,
1705 fe_degree,
1706 transpose>
1707 helper(given_degree,
1708 type_x,
1709 type_y,
1710 type_z,
1711 vv,
1712 interpolation_matrices,
1713 values);
1714
1715 if (faces > 0)
1716 switch (faces)
1717 {
1718 case 0:
1719 break;
1720 case 1:
1721 helper
1722 .template process_faces_fast<true, false, false>();
1723 break;
1724 case 2:
1725 helper
1726 .template process_faces_fast<false, true, false>();
1727 break;
1728 case 3:
1729 helper.template process_faces<true, true, false>();
1730 break;
1731 case 4:
1732 helper
1733 .template process_faces_fast<false, false, true>();
1734 break;
1735 case 5:
1736 helper.template process_faces<true, false, true>();
1737 break;
1738 case 6:
1739 helper.template process_faces<false, true, true>();
1740 break;
1741 case 7:
1742 helper.template process_faces<true, true, true>();
1743 break;
1744 }
1745
1746 if (edges > 0)
1747 switch (edges)
1748 {
1749 case 0:
1750 break;
1751 case 1:
1752 helper.template process_edge<true, false, false>();
1753 break;
1754 case 2:
1755 helper.template process_edge<false, true, false>();
1756 break;
1757 case 3:
1758 helper.template process_edge<true, true, false>();
1759 break;
1760 case 4:
1761 helper.template process_edge<false, false, true>();
1762 break;
1763 case 5:
1764 helper.template process_edge<true, false, true>();
1765 break;
1766 case 6:
1767 helper.template process_edge<false, true, true>();
1768 break;
1769 case 7:
1770 helper.template process_edge<true, true, true>();
1771 break;
1772 }
1773 }
1774 else
1775 {
1777 }
1778 }
1779
1780 values += shape_info.dofs_per_component_on_cell;
1781 }
1782 }
1783 };
1784
1785
1786
1787 template <int dim, typename Number>
1789 {
1790 public:
1791 template <int fe_degree, typename Number2>
1792 static bool
1793 run(const unsigned int n_desired_components,
1795 const bool transpose,
1797 Number::size()> &c_mask,
1798 Number *values)
1799 {
1800 using RunnerType =
1802 dim,
1803 fe_degree,
1804 Number>;
1805
1806 if (transpose)
1807 RunnerType::template run_internal<true>(n_desired_components,
1808 shape_info,
1809 c_mask,
1810 values);
1811 else
1812 RunnerType::template run_internal<false>(n_desired_components,
1813 shape_info,
1814 c_mask,
1815 values);
1816
1817 return false;
1818 }
1819
1820 template <int fe_degree>
1823 {
1824 return ((Number::size() > 2) && (fe_degree == -1 || fe_degree > 2)) ?
1827 }
1828 };
1829
1830
1831} // end of namespace internal
1832
1833#undef DEAL_II_ALWAYS_INLINE_RELEASE
1834
1835
1837
1838#endif
static void run_internal(const unsigned int n_desired_components, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const std::array< MatrixFreeFunctions::compressed_constraint_kind, Number::size()> &constraint_mask, Number *values)
static DEAL_II_ALWAYS_INLINE_RELEASE void interpolate_2D(const unsigned int given_degree, const typename Trait< Number, VectorizationType >::index_type v, const typename Trait< Number, VectorizationType >::interpolation_type *DEAL_II_RESTRICT weight, Number *DEAL_II_RESTRICT values)
static void run_internal(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const std::array< MatrixFreeFunctions::compressed_constraint_kind, Number::size()> &constraint_mask, Number *values)
static void interpolate(const unsigned int offset, const unsigned int outer_stride, const unsigned int given_degree, const Number mask_weight, const Number mask_write, const Number2 *DEAL_II_RESTRICT weights, Number *DEAL_II_RESTRICT values)
DEAL_II_ALWAYS_INLINE_RELEASE FEEvaluationImplHangingNodesScalarEntityInterpolationImpl(const unsigned int &given_degree, const bool &type_x, const bool &type_y, const bool &type_z, const typename Trait< Number, VectorizationType >::index_type &v, const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > &interpolation_matrices, Number *values)
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines_plane(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines_plane(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE FEEvaluationImplHangingNodesScalarEntityInterpolationImpl(const unsigned int &given_degree, const bool &type_x, const bool &type_y, const bool &type_z, const typename Trait< Number, VectorizationType >::index_type &v, const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > &interpolation_matrices, Number *values)
const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > & interpolation_matrices
static DEAL_II_ALWAYS_INLINE_RELEASE void interpolate_3D_face(const unsigned int dof_offset, const unsigned int given_degree, const typename Trait< Number, VectorizationType >::index_type v, const typename Trait< Number, VectorizationType >::interpolation_type *DEAL_II_RESTRICT weight, Number *DEAL_II_RESTRICT values)
DEAL_II_ALWAYS_INLINE_RELEASE FEEvaluationImplHangingNodesScalarEntityInterpolation(const T &t, const unsigned int &given_degree, const bool &type_x, const bool &type_y, const bool &type_z, const typename Trait< Number, VectorizationType >::index_type &v, const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > &interpolation_matrices, Number *values)
static DEAL_II_ALWAYS_INLINE_RELEASE void interpolate_3D_edge(const unsigned int p, const unsigned int given_degree, const typename Trait< Number, VectorizationType >::index_type v, const typename Trait< Number, VectorizationType >::interpolation_type *DEAL_II_RESTRICT weight, Number *DEAL_II_RESTRICT values)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_RESTRICT
Definition config.h:167
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition config.h:636
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition config.h:680
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_ALWAYS_INLINE_RELEASE
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
#define AssertIndexRange(index, range)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std::uint8_t compressed_constraint_kind
Definition dof_info.h:84
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:105
static bool run(const unsigned int n_desired_components, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, Number::size()> &c_mask, Number *values)
static constexpr FEEvaluationImplHangingNodesRunnerTypes used_runner_type()
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:490
static const std::array< AlignedVector< T1 >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE index_type create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 get_value(const T1 &value, const index_type &)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static const std::array< AlignedVector< interpolation_type >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE T1::value_type get_value(const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const typename T1::value_type &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE unsigned int create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE T1::value_type get_value(const typename T1::value_type &value, const index_type &i)
static const std::array< AlignedVector< T1 >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE index_type create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 get_value(const T1 &value, const index_type &)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static const std::array< AlignedVector< T1 >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 get_value(const T1 &value, const index_type &)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)