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