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