deal.II version GIT relicensing-2167-g9622207b8f 2024-11-21 12:40:00+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: LGPL-2.1-or-later
4// Copyright (C) 2021 - 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
16#ifndef dealii_matrix_free_evaluation_kernels_hanging_nodes_h
17#define dealii_matrix_free_evaluation_kernels_hanging_nodes_h
18
19#include <deal.II/base/config.h>
20
23
26
27
29
30#ifdef DEBUG
31# define DEAL_II_ALWAYS_INLINE_RELEASE
32#else
33# define DEAL_II_ALWAYS_INLINE_RELEASE DEAL_II_ALWAYS_INLINE
34#endif
35
36
37
38namespace internal
39{
45
46
47
53 {
57 index,
61 group,
66 mask,
70 sorted
71 };
72
73
74
76 int dim,
77 int fe_degree,
78 typename Number>
80
81
82
83 template <int dim, int fe_degree, typename Number>
86 dim,
87 fe_degree,
88 Number>
89 {
90 private:
91 template <int structdim,
92 unsigned int direction,
93 bool transpose,
94 typename Number2>
95 static void
96 interpolate(const unsigned int offset,
97 const unsigned int outer_stride,
98 const unsigned int given_degree,
99 const Number mask_weight,
100 const Number mask_write,
101 const Number2 *DEAL_II_RESTRICT weights,
102 Number *DEAL_II_RESTRICT values)
103 {
104 static constexpr unsigned int max_n_points_1D = 40;
105
106 static_assert(structdim == 1 || structdim == 2,
107 "Only 1D and 2d interpolation implemented");
108 Number temp[fe_degree != -1 ? fe_degree + 1 : max_n_points_1D];
109
110 const unsigned int points =
111 (fe_degree != -1 ? fe_degree : given_degree) + 1;
112
113 AssertIndexRange(points, max_n_points_1D);
114
115 const unsigned int stride = Utilities::pow(points, direction);
116
117 const unsigned int end_of_outer_loop = structdim == 1 ? 2 : points - 1;
118 for (unsigned int g = 1; g < end_of_outer_loop; ++g)
119 {
120 const unsigned int my_offset =
121 offset + (structdim > 1 ? g * outer_stride : 0);
122
123 // extract values to interpolate
124 for (unsigned int k = 0; k < points; ++k)
125 temp[k] = values[my_offset + k * stride];
126
127 // perform interpolation point by point and write back
128 for (unsigned int k = 0; k < points / 2; ++k)
129 {
130 const unsigned int kmirror = points - 1 - k;
131 Number sum0 = Number(), sum1 = Number(), sum2 = Number(),
132 sum3 = Number();
133 for (unsigned int h = 0; h < points; ++h)
134 {
135 const unsigned int hmirror = points - 1 - h;
136 // load from both sides of the interpolation matrix to
137 // reflect symmetry between the two subfaces along that
138 // direction
139 const Number w0 = weights[(transpose ? 1 : points) * kmirror +
140 (transpose ? points : 1) * hmirror];
141 const Number w1 = weights[(transpose ? 1 : points) * k +
142 (transpose ? points : 1) * h];
143 sum0 += temp[h] * w0;
144 sum1 += temp[h] * w1;
145 sum2 += temp[hmirror] * w1;
146 sum3 += temp[hmirror] * w0;
147 }
148 values[my_offset + k * stride] =
149 temp[k] +
150 mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
151 values[my_offset + kmirror * stride] =
152 temp[kmirror] +
153 mask_write *
154 (sum2 + mask_weight * (sum3 - sum2) - temp[kmirror]);
155 }
156
157 // cleanup case
158 if (points % 2)
159 {
160 const unsigned int k = points / 2;
161 Number sum0 = temp[k] * weights[(transpose ? 1 : points) * k +
162 (transpose ? points : 1) * k],
163 sum1 = sum0;
164 for (unsigned int h = 0; h < points / 2; ++h)
165 {
166 const unsigned int hmirror = points - 1 - h;
167 const Number w0 = weights[(transpose ? 1 : points) * k +
168 (transpose ? points : 1) * hmirror];
169 const Number w1 = weights[(transpose ? 1 : points) * k +
170 (transpose ? points : 1) * h];
171 sum0 += temp[h] * w0;
172 sum0 += temp[hmirror] * w1;
173 sum1 += temp[h] * w1;
174 sum1 += temp[hmirror] * w0;
175 }
176 values[my_offset + k * stride] =
177 temp[k] +
178 mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
179 }
180 }
181 }
182
183 public:
184 template <bool transpose, typename Number2>
185 static void
187 const unsigned int n_components,
190 Number::size()> &constraint_mask,
191 Number *values)
192 {
193 const unsigned int given_degree =
194 fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
195
196 const Number2 *DEAL_II_RESTRICT weights =
197 shape_info.data.front().subface_interpolation_matrices[0].data();
198
199 const unsigned int points = given_degree + 1;
200 const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
201
202 if (dim == 2)
203 {
204 ::ndarray<Number, 2> mask_weights = {};
205 ::ndarray<Number, 2, 2> mask_write = {};
206 ::ndarray<bool, 2, 2> do_face = {};
207
208 for (unsigned int v = 0; v < Number::size(); ++v)
209 {
210 const auto kind = constraint_mask[v];
211 const bool subcell_x = (kind >> 0) & 1;
212 const bool subcell_y = (kind >> 1) & 1;
213 const bool face_x = (kind >> 3) & 1;
214 const bool face_y = (kind >> 4) & 1;
215
216 if (face_y)
217 {
218 const unsigned int side = !subcell_y;
219 mask_write[0][side][v] = 1;
220 do_face[0][side] = true;
221 mask_weights[0][v] = subcell_x;
222 }
223
224 if (face_x)
225 {
226 const unsigned int side = !subcell_x;
227 mask_write[1][side][v] = 1;
228 do_face[1][side] = true;
229 mask_weights[1][v] = subcell_y;
230 }
231 }
232
233 // x direction
234 {
235 const std::array<unsigned int, 2> offsets = {
236 {0, (points - 1) * points}};
237 for (unsigned int c = 0; c < n_components; ++c)
238 for (unsigned int face = 0; face < 2; ++face)
239 if (do_face[0][face])
240 interpolate<1, 0, transpose>(offsets[face],
241 0,
242 given_degree,
243 mask_weights[0],
244 mask_write[0][face],
245 weights,
246 values + c * n_dofs);
247 }
248
249 // y direction
250 {
251 const std::array<unsigned int, 2> offsets = {{0, points - 1}};
252 for (unsigned int c = 0; c < n_components; ++c)
253 for (unsigned int face = 0; face < 2; ++face)
254 if (do_face[1][face])
255 interpolate<1, 1, transpose>(offsets[face],
256 0,
257 given_degree,
258 mask_weights[1],
259 mask_write[1][face],
260 weights,
261 values + c * n_dofs);
262 }
263 }
264 else if (dim == 3)
265 {
266 const unsigned int p0 = 0;
267 const unsigned int p1 = points - 1;
268 const unsigned int p2 = points * points - points;
269 const unsigned int p3 = points * points - 1;
270 const unsigned int p4 = points * points * points - points * points;
271 const unsigned int p5 =
272 points * points * points - points * points + points - 1;
273 const unsigned int p6 = points * points * points - points;
274
275 ::ndarray<bool, 3, 4> process_edge = {};
276 ::ndarray<bool, 3, 4> process_face = {};
277 ::ndarray<Number, 3, 4> mask_edge = {};
278 ::ndarray<Number, 3, 4> mask_face = {};
279 ::ndarray<Number, 3> mask_weights = {};
280
281 for (unsigned int v = 0; v < Number::size(); ++v)
282 {
283 const auto kind = constraint_mask[v];
284
285 const bool subcell_x = (kind >> 0) & 1;
286 const bool subcell_y = (kind >> 1) & 1;
287 const bool subcell_z = (kind >> 2) & 1;
288 const bool face_x = ((kind >> 3) & 1) ? (kind >> 5) & 1 : 0;
289 const bool face_y = ((kind >> 3) & 1) ? (kind >> 6) & 1 : 0;
290 const bool face_z = ((kind >> 3) & 1) ? (kind >> 7) & 1 : 0;
291 const bool edge_x = ((kind >> 4) & 1) ? (kind >> 5) & 1 : 0;
292 const bool edge_y = ((kind >> 4) & 1) ? (kind >> 6) & 1 : 0;
293 const bool edge_z = ((kind >> 4) & 1) ? (kind >> 7) & 1 : 0;
294
295 if (subcell_x)
296 mask_weights[0][v] = 1;
297 if (subcell_y)
298 mask_weights[1][v] = 1;
299 if (subcell_z)
300 mask_weights[2][v] = 1;
301
302 if (face_x)
303 {
304 const unsigned int side = !subcell_x;
305
306 mask_face[1][side][v] = process_face[1][side] = true;
307 mask_edge[1][side][v] = process_edge[1][side] = true;
308 mask_edge[1][2 + side][v] = process_edge[1][2 + side] = true;
309 mask_face[2][side][v] = process_face[2][side] = true;
310 mask_edge[2][side][v] = process_edge[2][side] = true;
311 mask_edge[2][2 + side][v] = process_edge[2][2 + side] = true;
312 }
313 if (face_y)
314 {
315 const unsigned int side = !subcell_y;
316
317 mask_face[0][side][v] = process_face[0][side] = true;
318 mask_edge[0][side][v] = process_edge[0][side] = true;
319 mask_edge[0][2 + side][v] = process_edge[0][2 + side] = true;
320 mask_face[2][2 + side][v] = process_face[2][2 + side] = true;
321 mask_edge[2][2 * side][v] = process_edge[2][2 * side] = true;
322 mask_edge[2][2 * side + 1][v] =
323 process_edge[2][2 * side + 1] = true;
324 }
325 if (face_z)
326 {
327 const unsigned int side = !subcell_z;
328
329 mask_face[0][2 + side][v] = process_face[0][2 + side] = true;
330 mask_edge[0][2 * side][v] = process_edge[0][2 * side] = true;
331 mask_edge[0][2 * side + 1][v] =
332 process_edge[0][2 * side + 1] = true;
333 mask_face[1][2 + side][v] = process_face[1][2 + side] = true;
334 mask_edge[1][2 * side][v] = process_edge[1][2 * side] = true;
335 mask_edge[1][2 * side + 1][v] =
336 process_edge[1][2 * side + 1] = true;
337 }
338 if (edge_x)
339 {
340 const unsigned int index = (!subcell_z) * 2 + (!subcell_y);
341 mask_edge[0][index][v] = process_edge[0][index] = true;
342 }
343 if (edge_y)
344 {
345 const unsigned int index = (!subcell_z) * 2 + (!subcell_x);
346 mask_edge[1][index][v] = process_edge[1][index] = true;
347 }
348 if (edge_z)
349 {
350 const unsigned int index = (!subcell_y) * 2 + (!subcell_x);
351 mask_edge[2][index][v] = process_edge[2][index] = true;
352 }
353 }
354
355 // direction 0:
356 if (given_degree > 1)
357 {
358 const std::array<unsigned int, 4> face_offsets = {
359 {p0, p2, p0, p4}};
360 const std::array<unsigned int, 2> outer_strides = {
361 {points * points, points}};
362 for (unsigned int c = 0; c < n_components; ++c)
363 for (unsigned int face = 0; face < 4; ++face)
364 if (process_face[0][face])
365 interpolate<2, 0, transpose>(face_offsets[face],
366 outer_strides[face / 2],
367 given_degree,
368 mask_weights[0],
369 mask_face[0][face],
370 weights,
371 values + c * n_dofs);
372 }
373 {
374 const std::array<unsigned int, 4> edge_offsets = {{p0, p2, p4, p6}};
375 for (unsigned int c = 0; c < n_components; ++c)
376 for (unsigned int edge = 0; edge < 4; ++edge)
377 if (process_edge[0][edge])
378 interpolate<1, 0, transpose>(edge_offsets[edge],
379 0,
380 given_degree,
381 mask_weights[0],
382 mask_edge[0][edge],
383 weights,
384 values + c * n_dofs);
385 }
386
387 // direction 1:
388 if (given_degree > 1)
389 {
390 const std::array<unsigned int, 4> face_offsets = {
391 {p0, p1, p0, p4}};
392 const std::array<unsigned int, 2> outer_strides = {
393 {points * points, 1}};
394 for (unsigned int c = 0; c < n_components; ++c)
395 for (unsigned int face = 0; face < 4; ++face)
396 if (process_face[1][face])
397 interpolate<2, 1, transpose>(face_offsets[face],
398 outer_strides[face / 2],
399 given_degree,
400 mask_weights[1],
401 mask_face[1][face],
402 weights,
403 values + c * n_dofs);
404 }
405
406 {
407 const std::array<unsigned int, 4> edge_offsets = {{p0, p1, p4, p5}};
408 for (unsigned int c = 0; c < n_components; ++c)
409 for (unsigned int edge = 0; edge < 4; ++edge)
410 if (process_edge[1][edge])
411 interpolate<1, 1, transpose>(edge_offsets[edge],
412 0,
413 given_degree,
414 mask_weights[1],
415 mask_edge[1][edge],
416 weights,
417 values + c * n_dofs);
418 }
419
420 // direction 2:
421 if (given_degree > 1)
422 {
423 const std::array<unsigned int, 4> face_offsets = {
424 {p0, p1, p0, p2}};
425 const std::array<unsigned int, 2> outer_strides = {{points, 1}};
426 for (unsigned int c = 0; c < n_components; ++c)
427 for (unsigned int face = 0; face < 4; ++face)
428 if (process_face[2][face])
429 interpolate<2, 2, transpose>(face_offsets[face],
430 outer_strides[face / 2],
431 given_degree,
432 mask_weights[2],
433 mask_face[2][face],
434 weights,
435 values + c * n_dofs);
436 }
437
438 {
439 const std::array<unsigned int, 4> edge_offsets = {{p0, p1, p2, p3}};
440 for (unsigned int c = 0; c < n_components; ++c)
441 for (unsigned int edge = 0; edge < 4; ++edge)
442 if (process_edge[2][edge])
443 interpolate<1, 2, transpose>(edge_offsets[edge],
444 0,
445 given_degree,
446 mask_weights[2],
447 mask_edge[2][edge],
448 weights,
449 values + c * n_dofs);
450 }
451 }
452 else
453 {
455 }
456 }
457 };
458
459 template <typename T1, VectorizationTypes VT>
460 struct Trait;
461
462 template <typename T1>
464 {
465 using value_type = typename T1::value_type;
466 using index_type = unsigned int;
468
469 template <typename T>
470 static inline const std::array<AlignedVector<interpolation_type>, 2> &
471 get_interpolation_matrix(const T &shape_info)
472 {
473 return shape_info.data.front().subface_interpolation_matrices_scalar;
474 }
475
476 static inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
478 T1::size()> mask,
480 T1::size()> mask_new,
481 const unsigned int v)
482 {
483 (void)mask;
484 (void)mask_new;
485 return v;
486 }
487
488 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
489 do_break(unsigned int v,
491 {
492 (void)v;
493 (void)kind;
494 return false;
495 }
496
497 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
498 do_continue(unsigned int v,
500 {
501 (void)v;
502 return kind ==
504 }
505
510 T1::size()> mask)
511 {
512 return mask;
513 }
514
515 static inline DEAL_II_ALWAYS_INLINE_RELEASE typename T1::value_type
516 get_value(const typename T1::value_type &value, const index_type &i)
517 {
518 (void)i;
519 return value;
520 }
521
522 static inline DEAL_II_ALWAYS_INLINE_RELEASE typename T1::value_type
523 get_value(const T1 &value, const index_type &i)
524 {
525 return value[i];
526 }
527
528 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
529 set_value(T1 &result,
530 const typename T1::value_type &value,
531 const index_type &i)
532 {
533 result[i] = value;
534 }
535 };
536
537 template <typename T1>
539 {
540 using value_type = T1;
541 using index_type = std::pair<T1, T1>;
543
544 template <typename T>
545 static inline const std::array<AlignedVector<T1>, 2> &
546 get_interpolation_matrix(const T &shape_info)
547 {
548 return shape_info.data.front().subface_interpolation_matrices;
549 }
550
551 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
552 do_break(unsigned int v,
554 {
555 (void)v;
556 (void)kind;
557 return false;
558 }
559
560 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
561 do_continue(unsigned int v,
563 {
564 (void)v;
565 return kind ==
567 }
568
569 static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type
571 T1::size()> mask,
573 T1::size()> mask_new,
574 const unsigned int v)
575 {
576 (void)mask;
577 (void)mask_new;
578 T1 result = 0.0;
579 result[v] = 1.0;
580 return {result, T1(1.0) - result};
581 }
582
587 T1::size()> mask)
588 {
589 return mask;
590 }
591
592 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
593 get_value(const T1 &value, const index_type &)
594 {
595 return value;
596 }
597
598 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
599 set_value(T1 &result, const T1 &value, const index_type &i)
600 {
601 result = result * i.second + value * i.first;
602 }
603 };
604
605 template <typename T1>
607 {
608 using value_type = T1;
609 using index_type = std::pair<T1, T1>;
611
612 template <typename T>
613 static inline const std::array<AlignedVector<T1>, 2> &
614 get_interpolation_matrix(const T &shape_info)
615 {
616 return shape_info.data.front().subface_interpolation_matrices;
617 }
618
619 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
620 do_break(unsigned int v,
622 {
623 (void)v;
624 return kind ==
626 }
627
628 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
629 do_continue(unsigned int v,
631 {
632 (void)v;
633 return kind ==
635 }
636
637 static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type
639 T1::size()> mask,
641 T1::size()> mask_new,
642 const unsigned int v)
643 {
644 T1 result;
645
646 for (unsigned int i = 0; i < T1::size(); ++i)
647 result[i] = mask_new[v] == mask[i];
648
649 return {result, T1(1.0) - result};
650 }
651
656 T1::size()> mask)
657 {
658 auto new_mask = mask;
659
660 std::sort(new_mask.begin(), new_mask.end());
661 std::fill(std::unique(new_mask.begin(), new_mask.end()),
662 new_mask.end(),
664
665 return new_mask;
666 }
667
668 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
669 get_value(const T1 &value, const index_type &)
670 {
671 return value;
672 }
673
674 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
675 set_value(T1 &result, const T1 &value, const index_type &i)
676 {
677 result = result * i.second + value * i.first;
678 }
679 };
680
681 template <typename T1>
683 {
684 using value_type = T1;
685 using index_type = T1;
687
688 template <typename T>
689 static inline const std::array<AlignedVector<T1>, 2> &
690 get_interpolation_matrix(const T &shape_info)
691 {
692 return shape_info.data.front().subface_interpolation_matrices;
693 }
694
695 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
696 do_break(unsigned int v,
698 {
699 (void)kind;
700 return v > 0;
701 }
702
703 static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
704 do_continue(unsigned int v,
706 {
707 (void)kind;
708
710
711 return v > 0; // should not be called
712 }
713
714 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
716 T1::size()> mask,
718 T1::size()> mask_new,
719 const unsigned int v)
720 {
721 (void)mask;
722 (void)mask_new;
723 (void)v;
724 return 1.0; // return something since not used
725 }
726
731 T1::size()> mask)
732 {
733 return mask;
734 }
735
736 static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
737 get_value(const T1 &value, const index_type &)
738 {
739 return value;
740 }
741
742 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
743 set_value(T1 &result, const T1 &value, const index_type &i)
744 {
745 (void)i;
746 result = value;
747 }
748 };
749
750
751
752 template <typename T,
753 typename Number,
754 VectorizationTypes VectorizationType,
755 int fe_degree,
756 bool transpose>
758 {
759 public:
762 const T &t,
763 const unsigned int &given_degree,
764 const bool &type_x,
765 const bool &type_y,
766 const bool &type_z,
768 const std::array<
772 Number *values)
773 : t(t)
775 , type_x(type_x)
776 , type_y(type_y)
777 , type_z(type_z)
778 , v(v)
780 , values(values)
781 {}
782
783 template <unsigned int direction, unsigned int d, bool skip_borders>
784 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
786 const unsigned int dof_offset,
787 const unsigned int given_degree,
790 *DEAL_II_RESTRICT weight,
791 Number *DEAL_II_RESTRICT values)
792 {
793 static constexpr unsigned int max_n_points_1D = 40;
794
796 temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
797
798 const unsigned int points =
799 (fe_degree != -1 ? fe_degree : given_degree) + 1;
800
801 AssertIndexRange(given_degree, max_n_points_1D);
802
803 const unsigned int stride = fe_degree != -1 ?
804 Utilities::pow(fe_degree + 1, direction) :
805 Utilities::pow(given_degree + 1, direction);
806
807 // direction side0 side1 side2
808 // 0 - p^2 p
809 // 1 p^2 - 1
810 // 2 p - 1
811 const unsigned int stride2 =
812 ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ?
813 (points * points) :
814 (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points :
815 1);
816
817 for (unsigned int g = (skip_borders ? 1 : 0);
818 g < points - (skip_borders ? 1 : 0);
819 ++g)
820 {
821 // copy result back
822 for (unsigned int k = 0; k < points; ++k)
824 values[dof_offset + k * stride + stride2 * g], v);
825
826 // perform interpolation point by point
827 for (unsigned int k = 0; k < points; ++k)
828 {
830 weight[(transpose ? 1 : points) * k], v) *
831 temp[0];
832 for (unsigned int h = 1; h < points; ++h)
834 weight[(transpose ? 1 : points) * k +
835 (transpose ? points : 1) * h],
836 v) *
837 temp[h];
839 values[dof_offset + k * stride + stride2 * g], sum, v);
840 }
841 }
842 }
843
844 template <unsigned int direction>
845 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
847 const unsigned int p,
848 const unsigned int given_degree,
851 *DEAL_II_RESTRICT weight,
852 Number *DEAL_II_RESTRICT values)
853 {
854 static constexpr unsigned int max_n_points_1D = 40;
855
857 temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
858
859 const unsigned int points =
860 (fe_degree != -1 ? fe_degree : given_degree) + 1;
861
862 AssertIndexRange(given_degree, max_n_points_1D);
863
864 const unsigned int stride = fe_degree != -1 ?
865 Utilities::pow(fe_degree + 1, direction) :
866 Utilities::pow(given_degree + 1, direction);
867
868 // copy result back
869 for (unsigned int k = 0; k < points; ++k)
870 temp[k] =
872 v);
873
874 // perform interpolation point by point
875 for (unsigned int k = 0; k < points; ++k)
876 {
878 weight[(transpose ? 1 : points) * k], v) *
879 temp[0];
880 for (unsigned int h = 1; h < points; ++h)
882 weight[(transpose ? 1 : points) * k +
883 (transpose ? points : 1) * h],
884 v) *
885 temp[h];
887 sum,
888 v);
889 }
890 }
891
892 template <bool do_x, bool do_y, bool do_z>
895 {
896 if (do_x)
897 interpolate_3D_edge<0>(t.line(0, type_y, type_z),
899 v,
901 values);
902
903 if (do_y)
904 interpolate_3D_edge<1>(t.line(1, type_x, type_z),
906 v,
908 values);
909
910 if (do_z)
911 interpolate_3D_edge<2>(t.line(2, type_x, type_y),
913 v,
915 values);
916 }
917
918 template <bool do_x, bool do_y, bool do_z>
921 {
922 static_assert((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) ||
923 (!do_x && !do_y && do_z),
924 "Only one face can be chosen.");
925
926 static const unsigned int direction = do_x ? 0 : (do_y ? 1 : 2);
927 const bool type = do_x ? type_x : (do_y ? type_y : type_z);
928
929 if (!do_x)
930 interpolate_3D_face<0, direction, false>(
931 t.face(direction, type),
933 v,
935 values);
936
937 if (!do_y)
938 interpolate_3D_face<1, direction, false>(
939 t.face(direction, type),
941 v,
943 values);
944
945 if (!do_z)
946 interpolate_3D_face<2, direction, false>(
947 t.face(direction, type),
949 v,
951 values);
952 }
953
954 template <bool do_x, bool do_y, bool do_z>
957 {
958 static_assert(((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) ||
959 (!do_x && !do_y && do_z)) == false,
960 "Only one face can be chosen.");
961
962 // direction 0
963 {
964 const auto inpterolation_matrix =
966
967 // faces
968 if (do_y && given_degree > 1)
969 interpolate_3D_face<0, 1, true>(
970 t.face(1, type_y), given_degree, v, inpterolation_matrix, values);
971
972 if (do_z && given_degree > 1)
973 interpolate_3D_face<0, 2, true>(
974 t.face(2, type_z), given_degree, v, inpterolation_matrix, values);
975
976 // direction 0 -> edges
977 interpolate_3D_edge<0>((do_x && do_y && !do_z) ?
978 (t.lines_plane(0, type_x, type_y, 0)) :
979 ((do_x && !do_y && do_z) ?
980 (t.lines_plane(1, type_x, type_z, 0)) :
981 (t.lines(0, type_y, type_z, 0))),
983 v,
984 inpterolation_matrix,
985 values);
986
987
988 interpolate_3D_edge<0>((do_x && do_y && !do_z) ?
989 (t.lines_plane(0, type_x, type_y, 1)) :
990 ((do_x && !do_y && do_z) ?
991 (t.lines_plane(1, type_x, type_z, 1)) :
992 (t.lines(0, type_y, type_z, 1))),
994 v,
995 inpterolation_matrix,
996 values);
997
998 if (do_y && do_z)
999 interpolate_3D_edge<0>(t.lines(0, type_y, type_z, 2),
1001 v,
1002 inpterolation_matrix,
1003 values);
1004 }
1005
1006 // direction 1
1007 {
1008 const auto inpterolation_matrix =
1010
1011 // faces
1012 if (do_x && given_degree > 1)
1013 interpolate_3D_face<1, 0, true>(
1014 t.face(0, type_x), given_degree, v, inpterolation_matrix, values);
1015
1016 if (do_z && given_degree > 1)
1017 interpolate_3D_face<1, 2, true>(
1018 t.face(2, type_z), given_degree, v, inpterolation_matrix, values);
1019
1020 // lines
1021 interpolate_3D_edge<1>((do_x && do_y && !do_z) ?
1022 (t.lines_plane(0, type_x, type_y, 2)) :
1023 ((!do_x && do_y && do_z) ?
1024 (t.lines_plane(2, type_y, type_z, 0)) :
1025 (t.lines(1, type_x, type_z, 0))),
1027 v,
1028 inpterolation_matrix,
1029 values);
1030
1031 interpolate_3D_edge<1>((do_x && do_y && !do_z) ?
1032 (t.lines_plane(0, type_x, type_y, 3)) :
1033 ((!do_x && do_y && do_z) ?
1034 (t.lines_plane(2, type_y, type_z, 1)) :
1035 (t.lines(1, type_x, type_z, 1))),
1037 v,
1038 inpterolation_matrix,
1039 values);
1040
1041 if (do_x && do_z)
1042 interpolate_3D_edge<1>(t.lines(1, type_x, type_z, 2),
1044 v,
1045 inpterolation_matrix,
1046 values);
1047 }
1048
1049 // direction 2 -> faces
1050 {
1051 const auto inpterolation_matrix =
1053
1054 if (do_x && given_degree > 1)
1055 interpolate_3D_face<2, 0, true>(
1056 t.face(0, type_x), given_degree, v, inpterolation_matrix, values);
1057
1058 if (do_y && given_degree > 1)
1059 interpolate_3D_face<2, 1, true>(
1060 t.face(1, type_y), given_degree, v, inpterolation_matrix, values);
1061
1062 // direction 2 -> edges
1063 interpolate_3D_edge<2>((do_x && !do_y && do_z) ?
1064 (t.lines_plane(1, type_x, type_z, 2)) :
1065 ((!do_x && do_y && do_z) ?
1066 (t.lines_plane(2, type_y, type_z, 2)) :
1067 (t.lines(2, type_x, type_y, 0))),
1069 v,
1070 inpterolation_matrix,
1071 values);
1072
1073 interpolate_3D_edge<2>((do_x && !do_y && do_z) ?
1074 (t.lines_plane(1, type_x, type_z, 3)) :
1075 ((!do_x && do_y && do_z) ?
1076 (t.lines_plane(2, type_y, type_z, 3)) :
1077 (t.lines(2, type_x, type_y, 1))),
1079 v,
1080 inpterolation_matrix,
1081 values);
1082
1083 if (do_x && do_y)
1084 interpolate_3D_edge<2>(t.lines(2, type_x, type_y, 2),
1086 v,
1087 inpterolation_matrix,
1088 values);
1089 }
1090 }
1091
1092 private:
1093 const T &t;
1094 const unsigned int &given_degree;
1095 const bool &type_x;
1096 const bool &type_y;
1097 const bool &type_z;
1099 const std::array<
1103 Number *values;
1104 };
1105
1109 enum class HelperType
1110 {
1115 constant,
1120 dynamic
1121 };
1122
1123 template <HelperType helper_type,
1124 typename Number,
1125 VectorizationTypes VectorizationType,
1126 int fe_degree,
1127 bool transpose>
1128 class Helper;
1129
1130 template <typename Number,
1131 VectorizationTypes VectorizationType,
1132 int fe_degree,
1133 bool transpose>
1135 Number,
1136 VectorizationType,
1137 fe_degree,
1138 transpose> : public HelperBase<Helper<HelperType::dynamic,
1139 Number,
1140 VectorizationType,
1141 fe_degree,
1142 transpose>,
1143 Number,
1144 VectorizationType,
1145 fe_degree,
1146 transpose>
1147 {
1148 public:
1150 Helper(const unsigned int &given_degree,
1151 const bool &type_x,
1152 const bool &type_y,
1153 const bool &type_z,
1155 const std::array<
1158 2> &interpolation_matrices,
1159 Number *values)
1161 Number,
1162 VectorizationType,
1163 fe_degree,
1164 transpose>,
1165 Number,
1166 VectorizationType,
1167 fe_degree,
1168 transpose>(*this,
1169 given_degree,
1170 type_x,
1171 type_y,
1172 type_z,
1173 v,
1174 interpolation_matrices,
1175 values)
1176 , points(given_degree + 1)
1177 {
1178 static_assert(fe_degree == -1, "Only working for fe_degree = -1.");
1179 }
1180
1181 const unsigned int points;
1182
1183 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1184 line(unsigned int i, unsigned int j, unsigned int k) const
1185 {
1186 return line_array[i][j][k];
1187 }
1188
1189 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1190 face(unsigned int i, unsigned int j) const
1191 {
1192 return face_array[i][j];
1193 }
1194
1195 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1196 lines_plane(unsigned int i,
1197 unsigned int j,
1198 unsigned int k,
1199 unsigned int l) const
1200 {
1201 return lines_plane_array[i][j][k][l];
1202 }
1203
1204 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1205 lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
1206 {
1207 return lines_array[i][j][k][l];
1208 }
1209
1210 private:
1211 const ::ndarray<unsigned int, 3, 2, 2> line_array = {
1212 {{{{{points * points * points - points, points *points - points}},
1213 {{points * points * points - points * points, 0}}}},
1214 {{{{points * points * points - points * points + points - 1,
1215 points - 1}},
1216 {{points * points * points - points * points, 0}}}},
1217 {{{{points * points - 1, points - 1}},
1218 {{points * points - points, 0}}}}}};
1219
1220 const ::ndarray<unsigned int, 3, 2> face_array = {
1221 {{{points - 1, 0}},
1222 {{points * points - points, 0}},
1223 {{points * points * points - points * points, 0}}}};
1224
1225 const ::ndarray<unsigned int, 3, 2, 2, 4> lines_plane_array = {
1226 {{{{{{{points * points - points,
1227 points *points *points - points,
1228 points - 1,
1229 points *points *points - points *points + points - 1}},
1230 {{0,
1231 points *points *points - points *points,
1232 points - 1,
1233 points *points *points - points *points + points - 1}}}},
1234 {{{{points * points - points,
1235 points *points *points - points,
1236 0,
1237 points *points *points - points *points}},
1238 {{0,
1239 points *points *points - points *points,
1240 0,
1241 points *points *points - points *points}}}}}},
1242 {{{{{{points * points * points - points * points,
1243 points *points *points - points,
1244 points - 1,
1245 points *points - 1}},
1246 {{0, points *points - points, points - 1, points *points - 1}}}},
1247 {{{{points * points * points - points * points,
1248 points *points *points - points,
1249 0,
1250 points *points - points}},
1251 {{0, points *points - points, 0, points *points - points}}}}}},
1252 {{{{{{points * points * points - points * points,
1253 points *points *points - points *points + points - 1,
1254 points *points - points,
1255 points *points - 1}},
1256 {{0, points - 1, points *points - points, points *points - 1}}}},
1257 {{{{points * points * points - points * points,
1258 points *points *points - points *points + points - 1,
1259 0,
1260 points - 1}},
1261 {{0, points - 1, 0, points - 1}}}}}}}};
1262
1263 const ::ndarray<unsigned int, 3, 2, 2, 3> lines_array = {
1264 {{{{{{{points * points - points,
1265 points *points *points - points *points,
1266 points *points *points - points}},
1267 {{0, points *points - points, points *points *points - points}}}},
1268 {{{{0,
1269 points *points *points - points *points,
1270 points *points *points - points}},
1271 {{0,
1272 points *points - points,
1273 points *points *points - points *points}}}}}},
1274 {{{{{{points - 1,
1275 points *points *points - points *points,
1276 points *points *points - points *points + points - 1}},
1277 {{0,
1278 points - 1,
1279 points *points *points - points *points + points - 1}}}},
1280 {{{{0,
1281 points *points *points - points *points,
1282 points *points *points - points *points + points - 1}},
1283 {{0, points - 1, points *points *points - points *points}}}}}},
1284 {{{{{{points - 1, points *points - points, points *points - 1}},
1285 {{0, points - 1, points *points - 1}}}},
1286 {{{{0, points *points - points, points *points - 1}},
1287 {{0, points - 1, points *points - points}}}}}}}};
1288 };
1289
1290 template <typename Number,
1291 VectorizationTypes VectorizationType,
1292 int fe_degree,
1293 bool transpose>
1295 Number,
1296 VectorizationType,
1297 fe_degree,
1298 transpose> : public HelperBase<Helper<HelperType::constant,
1299 Number,
1300 VectorizationType,
1301 fe_degree,
1302 transpose>,
1303 Number,
1304 VectorizationType,
1305 fe_degree,
1306 transpose>
1307 {
1308 public:
1310 Helper(const unsigned int &given_degree,
1311 const bool &type_x,
1312 const bool &type_y,
1313 const bool &type_z,
1315 const std::array<
1318 2> &interpolation_matrices,
1319 Number *values)
1321 Number,
1322 VectorizationType,
1323 fe_degree,
1324 transpose>,
1325 Number,
1326 VectorizationType,
1327 fe_degree,
1328 transpose>(*this,
1329 given_degree,
1330 type_x,
1331 type_y,
1332 type_z,
1333 v,
1334 interpolation_matrices,
1335 values)
1336 {
1337 static_assert(fe_degree != -1, "Only working for fe_degree != -1.");
1338 }
1339
1340
1341 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1342 line(unsigned int i, unsigned int j, unsigned int k) const
1343 {
1344 static constexpr unsigned int points = fe_degree + 1;
1345
1346 static constexpr ::ndarray<unsigned int, 3, 2, 2> line_array = {
1347 {{{{{points * points * points - points, points * points - points}},
1348 {{points * points * points - points * points, 0}}}},
1349 {{{{points * points * points - points * points + points - 1,
1350 points - 1}},
1351 {{points * points * points - points * points, 0}}}},
1352 {{{{points * points - 1, points - 1}},
1353 {{points * points - points, 0}}}}}};
1354
1355 return line_array[i][j][k];
1356 }
1357
1358 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1359 face(unsigned int i, unsigned int j) const
1360 {
1361 static constexpr unsigned int points = fe_degree + 1;
1362
1363 static constexpr ::ndarray<unsigned int, 3, 2> face_array = {
1364 {{{points - 1, 0}},
1365 {{points * points - points, 0}},
1366 {{points * points * points - points * points, 0}}}};
1367
1368 return face_array[i][j];
1369 }
1370
1371 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1372 lines_plane(unsigned int i,
1373 unsigned int j,
1374 unsigned int k,
1375 unsigned int l) const
1376 {
1377 static constexpr unsigned int points = fe_degree + 1;
1378
1379 static constexpr ::ndarray<unsigned int, 3, 2, 2, 4>
1380 lines_plane_array = {
1381 {{{{{{{points * points - points,
1382 points * points * points - points,
1383 points - 1,
1384 points * points * points - points * points + points - 1}},
1385 {{0,
1386 points * points * points - points * points,
1387 points - 1,
1388 points * points * points - points * points + points - 1}}}},
1389 {{{{points * points - points,
1390 points * points * points - points,
1391 0,
1392 points * points * points - points * points}},
1393 {{0,
1394 points * points * points - points * points,
1395 0,
1396 points * points * points - points * points}}}}}},
1397 {{{{{{points * points * points - points * points,
1398 points * points * points - points,
1399 points - 1,
1400 points * points - 1}},
1401 {{0,
1402 points * points - points,
1403 points - 1,
1404 points * points - 1}}}},
1405 {{{{points * points * points - points * points,
1406 points * points * points - points,
1407 0,
1408 points * points - points}},
1409 {{0, points * points - points, 0, points * points - points}}}}}},
1410 {{{{{{points * points * points - points * points,
1411 points * points * points - points * points + points - 1,
1412 points * points - points,
1413 points * points - 1}},
1414 {{0,
1415 points - 1,
1416 points * points - points,
1417 points * points - 1}}}},
1418 {{{{points * points * points - points * points,
1419 points * points * points - points * points + points - 1,
1420 0,
1421 points - 1}},
1422 {{0, points - 1, 0, points - 1}}}}}}}};
1423
1424 return lines_plane_array[i][j][k][l];
1425 }
1426
1427 inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1428 lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
1429 {
1430 static constexpr unsigned int points = fe_degree + 1;
1431
1432 static constexpr ::ndarray<unsigned int, 3, 2, 2, 3> lines_array = {
1433 {{{{{{{points * points - points,
1434 points * points * points - points * points,
1435 points * points * points - points}},
1436 {{0,
1437 points * points - points,
1438 points * points * points - points}}}},
1439 {{{{0,
1440 points * points * points - points * points,
1441 points * points * points - points}},
1442 {{0,
1443 points * points - points,
1444 points * points * points - points * points}}}}}},
1445 {{{{{{points - 1,
1446 points * points * points - points * points,
1447 points * points * points - points * points + points - 1}},
1448 {{0,
1449 points - 1,
1450 points * points * points - points * points + points - 1}}}},
1451 {{{{0,
1452 points * points * points - points * points,
1453 points * points * points - points * points + points - 1}},
1454 {{0, points - 1, points * points * points - points * points}}}}}},
1455 {{{{{{points - 1, points * points - points, points * points - 1}},
1456 {{0, points - 1, points * points - 1}}}},
1457 {{{{0, points * points - points, points * points - 1}},
1458 {{0, points - 1, points * points - points}}}}}}}};
1459
1460 return lines_array[i][j][k][l];
1461 }
1462 };
1463
1464
1465 template <int dim, int fe_degree, typename Number>
1468 dim,
1469 fe_degree,
1470 Number>
1471 {
1472 public:
1473 static const VectorizationTypes VectorizationType =
1475
1476 private:
1477 template <unsigned int side, bool transpose>
1478 static inline DEAL_II_ALWAYS_INLINE_RELEASE void
1480 const unsigned int given_degree,
1483 *DEAL_II_RESTRICT weight,
1484 Number *DEAL_II_RESTRICT values)
1485 {
1486 static constexpr unsigned int max_n_points_1D = 40;
1487
1489 temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
1490
1491 const unsigned int points =
1492 (fe_degree != -1 ? fe_degree : given_degree) + 1;
1493
1494 AssertIndexRange(given_degree, max_n_points_1D);
1495
1496 const unsigned int d = side / 2; // direction
1497 const unsigned int s = side % 2; // left or right surface
1498
1499 const unsigned int offset = ::Utilities::pow(points, d + 1);
1500 const unsigned int stride =
1501 (s == 0 ? 0 : (points - 1)) * ::Utilities::pow(points, d);
1502
1503 const unsigned int r1 = ::Utilities::pow(points, dim - d - 1);
1504 const unsigned int r2 = ::Utilities::pow(points, d);
1505
1506 // copy result back
1507 for (unsigned int i = 0, k = 0; i < r1; ++i)
1508 for (unsigned int j = 0; j < r2; ++j, ++k)
1510 values[i * offset + stride + j], v);
1511
1512 // perform interpolation point by point (note: r1 * r2 ==
1513 // points^(dim-1))
1514 for (unsigned int i = 0, k = 0; i < r1; ++i)
1515 for (unsigned int j = 0; j < r2; ++j, ++k)
1516 {
1518 for (unsigned int h = 0; h < points; ++h)
1520 weight[(transpose ? 1 : points) * k +
1521 (transpose ? points : 1) * h],
1522 v) *
1523 temp[h];
1525 values[i * offset + stride + j], sum, v);
1526 }
1527 }
1528
1529 public:
1530 template <bool transpose, typename Number2>
1531 static void
1533 const unsigned int n_desired_components,
1536 Number::size()> &constraint_mask,
1537 Number *values)
1538 {
1539 const unsigned int given_degree =
1540 fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
1541
1542 const auto &interpolation_matrices =
1544
1545 const auto constraint_mask_sorted =
1547
1548 for (unsigned int c = 0; c < n_desired_components; ++c)
1549 {
1550 for (unsigned int v = 0; v < Number::size(); ++v)
1551 {
1552 const auto mask = constraint_mask_sorted[v];
1553
1555 break;
1556
1558 continue;
1559
1560 const auto vv =
1562 constraint_mask_sorted,
1563 v);
1564
1565 if (dim == 2) // 2d: only faces
1566 {
1567 const bool subcell_x = (mask >> 0) & 1;
1568 const bool subcell_y = (mask >> 1) & 1;
1569 const bool face_x = (mask >> 3) & 1;
1570 const bool face_y = (mask >> 4) & 1;
1571
1572 // direction 0:
1573 if (face_y)
1574 {
1575 const auto *weights =
1576 interpolation_matrices[!subcell_x].data();
1577
1578 if (subcell_y)
1579 interpolate_2D<2, transpose>(given_degree,
1580 vv,
1581 weights,
1582 values); // face 2
1583 else
1584 interpolate_2D<3, transpose>(given_degree,
1585 vv,
1586 weights,
1587 values); // face 3
1588 }
1589
1590 // direction 1:
1591 if (face_x)
1592 {
1593 const auto *weights =
1594 interpolation_matrices[!subcell_y].data();
1595
1596 if (subcell_x)
1597 interpolate_2D<0, transpose>(given_degree,
1598 vv,
1599 weights,
1600 values); // face 0
1601 else
1602 interpolate_2D<1, transpose>(given_degree,
1603 vv,
1604 weights,
1605 values); // face 1
1606 }
1607 }
1608 else if (dim == 3) // 3d faces and edges
1609 {
1610 const bool type_x = (mask >> 0) & 1;
1611 const bool type_y = (mask >> 1) & 1;
1612 const bool type_z = (mask >> 2) & 1;
1613
1614 const auto flag_0 = (mask >> 3) & 3;
1615 const auto flag_1 = (mask >> 5) & 7;
1616 const auto faces = (flag_0 & 0b01) ? flag_1 : 0;
1617 const auto edges = (flag_0 & 0b10) ? flag_1 : 0;
1618
1619 Helper<fe_degree == -1 ? HelperType::dynamic :
1621 Number,
1622 VectorizationType,
1623 fe_degree,
1624 transpose>
1625 helper(given_degree,
1626 type_x,
1627 type_y,
1628 type_z,
1629 vv,
1630 interpolation_matrices,
1631 values);
1632
1633 if (faces > 0)
1634 switch (faces)
1635 {
1636 case 0:
1637 break;
1638 case 1:
1639 helper
1640 .template process_faces_fast<true, false, false>();
1641 break;
1642 case 2:
1643 helper
1644 .template process_faces_fast<false, true, false>();
1645 break;
1646 case 3:
1647 helper.template process_faces<true, true, false>();
1648 break;
1649 case 4:
1650 helper
1651 .template process_faces_fast<false, false, true>();
1652 break;
1653 case 5:
1654 helper.template process_faces<true, false, true>();
1655 break;
1656 case 6:
1657 helper.template process_faces<false, true, true>();
1658 break;
1659 case 7:
1660 helper.template process_faces<true, true, true>();
1661 break;
1662 }
1663
1664 if (edges > 0)
1665 switch (edges)
1666 {
1667 case 0:
1668 break;
1669 case 1:
1670 helper.template process_edge<true, false, false>();
1671 break;
1672 case 2:
1673 helper.template process_edge<false, true, false>();
1674 break;
1675 case 3:
1676 helper.template process_edge<true, true, false>();
1677 break;
1678 case 4:
1679 helper.template process_edge<false, false, true>();
1680 break;
1681 case 5:
1682 helper.template process_edge<true, false, true>();
1683 break;
1684 case 6:
1685 helper.template process_edge<false, true, true>();
1686 break;
1687 case 7:
1688 helper.template process_edge<true, true, true>();
1689 break;
1690 }
1691 }
1692 else
1693 {
1695 }
1696 }
1697
1698 values += shape_info.dofs_per_component_on_cell;
1699 }
1700 }
1701 };
1702
1703
1704
1705 template <int dim, typename Number>
1707 {
1708 public:
1709 template <int fe_degree, typename Number2>
1710 static bool
1711 run(const unsigned int n_desired_components,
1713 const bool transpose,
1715 Number::size()> &c_mask,
1716 Number *values)
1717 {
1718 using RunnerType =
1720 dim,
1721 fe_degree,
1722 Number>;
1723
1724 if (transpose)
1725 RunnerType::template run_internal<true>(n_desired_components,
1726 shape_info,
1727 c_mask,
1728 values);
1729 else
1730 RunnerType::template run_internal<false>(n_desired_components,
1731 shape_info,
1732 c_mask,
1733 values);
1734
1735 return false;
1736 }
1737
1738 template <int fe_degree>
1741 {
1742 return ((Number::size() > 2) && (fe_degree == -1 || fe_degree > 2)) ?
1745 }
1746 };
1747
1748
1749} // end of namespace internal
1750
1751#undef DEAL_II_ALWAYS_INLINE_RELEASE
1752
1753
1755
1756#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 void process_faces() const
const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > & interpolation_matrices
DEAL_II_ALWAYS_INLINE_RELEASE HelperBase(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)
const Trait< Number, VectorizationType >::index_type & v
DEAL_II_ALWAYS_INLINE_RELEASE void process_edge() const
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 void process_faces_fast() 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 Helper(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 face(unsigned int i, unsigned int j) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int line(unsigned int i, unsigned int j, unsigned int k) const
DEAL_II_ALWAYS_INLINE_RELEASE Helper(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_plane(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int face(unsigned int i, unsigned int j) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int line(unsigned int i, unsigned int j, unsigned int k) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_RESTRICT
Definition config.h:110
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_ALWAYS_INLINE_RELEASE
#define AssertIndexRange(index, range)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std::uint8_t compressed_constraint_kind
Definition dof_info.h:86
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
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:485
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)