deal.II version GIT relicensing-2287-g6548a49e0a 2024-12-20 18:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
evaluation_kernels_face.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 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_face_h
17#define dealii_matrix_free_evaluation_kernels_face_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/geometry_info.h>
22#include <deal.II/base/ndarray.h>
23#include <deal.II/base/utilities.h>
24#include <deal.II/base/vectorization.h>
25
26#include <deal.II/matrix_free/dof_info.h>
27#include <deal.II/matrix_free/evaluation_flags.h>
28#include <deal.II/matrix_free/fe_evaluation_data.h>
29#include <deal.II/matrix_free/shape_info.h>
30#include <deal.II/matrix_free/tensor_product_kernels.h>
31
32
33DEAL_II_NAMESPACE_OPEN
34
35
36namespace internal
37{
38 template <bool symmetric_evaluate,
39 int dim,
40 int fe_degree,
41 int n_q_points_1d,
42 typename Number>
43 struct FEFaceEvaluationImpl
44 {
45 // We enable a transformation to collocation for derivatives if it gives
46 // correct results (first two conditions), if it is the most efficient
47 // choice in terms of operation counts (third condition) and if we were
48 // able to initialize the fields in shape_info.templates.h from the
49 // polynomials (fourth condition).
50 using Number2 =
51 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
52
53 using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
54 evaluate_general,
55 dim - 1,
56 fe_degree + 1,
57 n_q_points_1d,
58 Number,
59 Number2>;
60
61 static Eval
62 create_evaluator_tensor_product(
63 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
64 const unsigned int subface_index,
65 const unsigned int direction)
66 {
67 if (symmetric_evaluate)
68 return Eval(data.shape_values_eo,
69 data.shape_gradients_eo,
70 data.shape_hessians_eo,
71 data.fe_degree + 1,
72 data.n_q_points_1d);
73 else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
74 return Eval(data.shape_values,
75 data.shape_gradients,
76 data.shape_hessians,
77 data.fe_degree + 1,
78 data.n_q_points_1d);
79 else
80 {
81 const unsigned int index =
82 direction == 0 ? subface_index % 2 : subface_index / 2;
83 return Eval(data.values_within_subface[index],
84 data.gradients_within_subface[index],
85 data.hessians_within_subface[index],
86 data.fe_degree + 1,
87 data.n_q_points_1d);
88 }
89 }
90
91 static void
92 evaluate_in_face(
93 const unsigned int n_components,
94 const EvaluationFlags::EvaluationFlags evaluation_flag,
95 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
96 Number *values_dofs,
97 Number *values_quad,
98 Number *gradients_quad,
99 Number *hessians_quad,
100 Number *scratch_data,
101 const unsigned int subface_index)
102 {
103 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
104 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
105
106 const std::size_t n_dofs = fe_degree > -1 ?
107 Utilities::pow(fe_degree + 1, dim - 1) :
108 Utilities::pow(data.fe_degree + 1, dim - 1);
109 const std::size_t n_q_points =
110 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
111 Utilities::pow(data.n_q_points_1d, dim - 1);
112
113 // keep a copy of the original pointer for the case of the Hessians
114 Number *values_dofs_ptr = values_dofs;
115
116 if ((evaluation_flag & EvaluationFlags::values) != 0u &&
117 ((evaluation_flag & EvaluationFlags::gradients) == 0u))
118 for (unsigned int c = 0; c < n_components; ++c)
119 {
120 switch (dim)
121 {
122 case 3:
123 eval0.template values<0, true, false>(values_dofs,
124 values_quad);
125 eval1.template values<1, true, false>(values_quad,
126 values_quad);
127 break;
128 case 2:
129 eval0.template values<0, true, false>(values_dofs,
130 values_quad);
131 break;
132 case 1:
133 values_quad[0] = values_dofs[0];
134 break;
135 default:
136 DEAL_II_NOT_IMPLEMENTED();
137 }
138 // Note: we always keep storage of values, 1st and 2nd derivatives
139 // in an array
140 values_dofs += 3 * n_dofs;
141 values_quad += n_q_points;
142 }
143 else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
144 for (unsigned int c = 0; c < n_components; ++c)
145 {
146 switch (dim)
147 {
148 case 3:
149 if (symmetric_evaluate &&
150 use_collocation_evaluation(fe_degree, n_q_points_1d))
151 {
152 eval0.template values<0, true, false>(values_dofs,
153 values_quad);
154 eval0.template values<1, true, false>(values_quad,
155 values_quad);
156 EvaluatorTensorProduct<evaluate_evenodd,
157 dim - 1,
158 n_q_points_1d,
159 n_q_points_1d,
160 Number,
161 Number2>
162 eval_grad({}, data.shape_gradients_collocation_eo, {});
163 eval_grad.template gradients<0, true, false, 3>(
164 values_quad, gradients_quad);
165 eval_grad.template gradients<1, true, false, 3>(
166 values_quad, gradients_quad + 1);
167 }
168 else
169 {
170 // grad x
171 eval0.template gradients<0, true, false>(values_dofs,
172 scratch_data);
173 eval1.template values<1, true, false, 3>(scratch_data,
174 gradients_quad);
175
176 // grad y
177 eval0.template values<0, true, false>(values_dofs,
178 scratch_data);
179 eval1.template gradients<1, true, false, 3>(
180 scratch_data, gradients_quad + 1);
181
182 if ((evaluation_flag & EvaluationFlags::values) != 0u)
183 eval1.template values<1, true, false>(scratch_data,
184 values_quad);
185 }
186 // grad z
187 eval0.template values<0, true, false>(values_dofs + n_dofs,
188 scratch_data);
189 eval1.template values<1, true, false, 3>(scratch_data,
190 gradients_quad + 2);
191
192 break;
193 case 2:
194 eval0.template values<0, true, false, 2>(values_dofs + n_dofs,
195 gradients_quad + 1);
196 eval0.template gradients<0, true, false, 2>(values_dofs,
197 gradients_quad);
198 if ((evaluation_flag & EvaluationFlags::values) != 0u)
199 eval0.template values<0, true, false>(values_dofs,
200 values_quad);
201 break;
202 case 1:
203 values_quad[0] = values_dofs[0];
204 gradients_quad[0] = values_dofs[1];
205 break;
206 default:
207 AssertThrow(false, ExcNotImplemented());
208 }
209 values_dofs += 3 * n_dofs;
210 values_quad += n_q_points;
211 gradients_quad += dim * n_q_points;
212 }
213
214 if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
215 {
216 values_dofs = values_dofs_ptr;
217 for (unsigned int c = 0; c < n_components; ++c)
218 {
219 switch (dim)
220 {
221 case 3:
222 // grad xx
223 eval0.template hessians<0, true, false>(values_dofs,
224 scratch_data);
225 eval1.template values<1, true, false>(scratch_data,
226 hessians_quad);
227
228 // grad yy
229 eval0.template values<0, true, false>(values_dofs,
230 scratch_data);
231 eval1.template hessians<1, true, false>(scratch_data,
232 hessians_quad +
233 n_q_points);
234
235 // grad zz
236 eval0.template values<0, true, false>(values_dofs +
237 2 * n_dofs,
238 scratch_data);
239 eval1.template values<1, true, false>(scratch_data,
240 hessians_quad +
241 2 * n_q_points);
242
243 // grad xy
244 eval0.template gradients<0, true, false>(values_dofs,
245 scratch_data);
246 eval1.template gradients<1, true, false>(scratch_data,
247 hessians_quad +
248 3 * n_q_points);
249
250 // grad xz
251 eval0.template gradients<0, true, false>(values_dofs +
252 n_dofs,
253 scratch_data);
254 eval1.template values<1, true, false>(scratch_data,
255 hessians_quad +
256 4 * n_q_points);
257
258 // grad yz
259 eval0.template values<0, true, false>(values_dofs + n_dofs,
260 scratch_data);
261 eval1.template gradients<1, true, false>(scratch_data,
262 hessians_quad +
263 5 * n_q_points);
264
265 break;
266 case 2:
267 // grad xx
268 eval0.template hessians<0, true, false>(values_dofs,
269 hessians_quad);
270 // grad yy
271 eval0.template values<0, true, false>(
272 values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
273 // grad xy
274 eval0.template gradients<0, true, false>(
275 values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
276 break;
277 case 1:
278 hessians_quad[0] = values_dofs[2];
279 break;
280 default:
281 AssertThrow(false, ExcNotImplemented());
282 }
283 values_dofs += 3 * n_dofs;
284 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
285 }
286 }
287 }
288
289 static void
290 integrate_in_face(
291 const unsigned int n_components,
292 const EvaluationFlags::EvaluationFlags integration_flag,
293 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
294 Number *values_dofs,
295 Number *values_quad,
296 Number *gradients_quad,
297 Number *hessians_quad,
298 Number *scratch_data,
299 const unsigned int subface_index)
300 {
301 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
302 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
303
304 const std::size_t n_dofs =
305 fe_degree > -1 ?
306 Utilities::pow(fe_degree + 1, dim - 1) :
307 (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
308 const std::size_t n_q_points =
309 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
310 Utilities::pow(data.n_q_points_1d, dim - 1);
311
312 // keep a copy of the original pointer for the case of the Hessians
313 Number *values_dofs_ptr = values_dofs;
314
315 if ((integration_flag & EvaluationFlags::values) != 0u &&
316 (integration_flag & EvaluationFlags::gradients) == 0u)
317 for (unsigned int c = 0; c < n_components; ++c)
318 {
319 switch (dim)
320 {
321 case 3:
322 eval1.template values<1, false, false>(values_quad,
323 values_quad);
324 eval0.template values<0, false, false>(values_quad,
325 values_dofs);
326 break;
327 case 2:
328 eval0.template values<0, false, false>(values_quad,
329 values_dofs);
330 break;
331 case 1:
332 values_dofs[0] = values_quad[0];
333 break;
334 default:
335 DEAL_II_NOT_IMPLEMENTED();
336 }
337 values_dofs += 3 * n_dofs;
338 values_quad += n_q_points;
339 }
340 else if ((integration_flag & EvaluationFlags::gradients) != 0u)
341 for (unsigned int c = 0; c < n_components; ++c)
342 {
343 switch (dim)
344 {
345 case 3:
346 // grad z
347 eval1.template values<1, false, false, 3>(gradients_quad + 2,
348 scratch_data);
349 eval0.template values<0, false, false>(scratch_data,
350 values_dofs + n_dofs);
351 if (symmetric_evaluate &&
352 use_collocation_evaluation(fe_degree, n_q_points_1d))
353 {
354 EvaluatorTensorProduct<evaluate_evenodd,
355 dim - 1,
356 n_q_points_1d,
357 n_q_points_1d,
358 Number,
359 Number2>
360 eval_grad({}, data.shape_gradients_collocation_eo, {});
361 if ((integration_flag & EvaluationFlags::values) != 0u)
362 eval_grad.template gradients<1, false, true, 3>(
363 gradients_quad + 1, values_quad);
364 else
365 eval_grad.template gradients<1, false, false, 3>(
366 gradients_quad + 1, values_quad);
367 eval_grad.template gradients<0, false, true, 3>(
368 gradients_quad, values_quad);
369 eval0.template values<1, false, false>(values_quad,
370 values_quad);
371 eval0.template values<0, false, false>(values_quad,
372 values_dofs);
373 }
374 else
375 {
376 if ((integration_flag & EvaluationFlags::values) != 0u)
377 {
378 eval1.template values<1, false, false>(values_quad,
379 scratch_data);
380 eval1.template gradients<1, false, true, 3>(
381 gradients_quad + 1, scratch_data);
382 }
383 else
384 eval1.template gradients<1, false, false, 3>(
385 gradients_quad + 1, scratch_data);
386
387 // grad y
388 eval0.template values<0, false, false>(scratch_data,
389 values_dofs);
390
391 // grad x
392 eval1.template values<1, false, false, 3>(gradients_quad,
393 scratch_data);
394 eval0.template gradients<0, false, true>(scratch_data,
395 values_dofs);
396 }
397 break;
398 case 2:
399 eval0.template values<0, false, false, 2>(gradients_quad + 1,
400 values_dofs +
401 n_dofs);
402 eval0.template gradients<0, false, false, 2>(gradients_quad,
403 values_dofs);
404 if ((integration_flag & EvaluationFlags::values) != 0u)
405 eval0.template values<0, false, true>(values_quad,
406 values_dofs);
407 break;
408 case 1:
409 values_dofs[0] = values_quad[0];
410 values_dofs[1] = gradients_quad[0];
411 break;
412 default:
413 AssertThrow(false, ExcNotImplemented());
414 }
415 values_dofs += 3 * n_dofs;
416 values_quad += n_q_points;
417 gradients_quad += dim * n_q_points;
418 }
419
420 if ((integration_flag & EvaluationFlags::hessians) != 0u)
421 {
422 values_dofs = values_dofs_ptr;
423 for (unsigned int c = 0; c < n_components; ++c)
424 {
425 switch (dim)
426 {
427 case 3:
428 // grad xx
429 eval1.template values<1, false, false>(hessians_quad,
430 scratch_data);
431 if ((integration_flag & (EvaluationFlags::values |
432 EvaluationFlags::gradients)) != 0u)
433 eval0.template hessians<0, false, true>(scratch_data,
434 values_dofs);
435 else
436 eval0.template hessians<0, false, false>(scratch_data,
437 values_dofs);
438
439 // grad yy
440 eval1.template hessians<1, false, false>(hessians_quad +
441 n_q_points,
442 scratch_data);
443 eval0.template values<0, false, true>(scratch_data,
444 values_dofs);
445
446 // grad zz
447 eval1.template values<1, false, false>(hessians_quad +
448 2 * n_q_points,
449 scratch_data);
450 eval0.template values<0, false, false>(scratch_data,
451 values_dofs +
452 2 * n_dofs);
453
454 // grad xy
455 eval1.template gradients<1, false, false>(hessians_quad +
456 3 * n_q_points,
457 scratch_data);
458 eval0.template gradients<0, false, true>(scratch_data,
459 values_dofs);
460
461 // grad xz
462 eval1.template values<1, false, false>(hessians_quad +
463 4 * n_q_points,
464 scratch_data);
465 if ((integration_flag & EvaluationFlags::gradients) != 0u)
466 eval0.template gradients<0, false, true>(scratch_data,
467 values_dofs +
468 n_dofs);
469 else
470 eval0.template gradients<0, false, false>(scratch_data,
471 values_dofs +
472 n_dofs);
473
474 // grad yz
475 eval1.template gradients<1, false, false>(hessians_quad +
476 5 * n_q_points,
477 scratch_data);
478 eval0.template values<0, false, true>(scratch_data,
479 values_dofs + n_dofs);
480
481 break;
482 case 2:
483 // grad xx
484 if ((integration_flag & (EvaluationFlags::values |
485 EvaluationFlags::gradients)) != 0u)
486 eval0.template hessians<0, false, true>(hessians_quad,
487 values_dofs);
488 else
489 eval0.template hessians<0, false, false>(hessians_quad,
490 values_dofs);
491
492 // grad yy
493 eval0.template values<0, false, false>(
494 hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
495 // grad xy
496 if ((integration_flag & EvaluationFlags::gradients) != 0u)
497 eval0.template gradients<0, false, true>(
498 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
499 else
500 eval0.template gradients<0, false, false>(
501 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
502 break;
503 case 1:
504 values_dofs[2] = hessians_quad[0];
505 if ((integration_flag & EvaluationFlags::values) == 0u)
506 values_dofs[0] = 0;
507 if ((integration_flag & EvaluationFlags::gradients) == 0u)
508 values_dofs[1] = 0;
509 break;
510 default:
511 AssertThrow(false, ExcNotImplemented());
512 }
513 values_dofs += 3 * n_dofs;
514 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
515 }
516 }
517 }
518 };
519
520
521
522 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
523 struct FEFaceEvaluationImplRaviartThomas
524 {
525 using Number2 =
526 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
527
532 template <bool do_integrate>
533 static inline void
534 evaluate_or_integrate_in_face(
535 const EvaluationFlags::EvaluationFlags evaluation_flag,
536 const std::vector<MatrixFreeFunctions::UnivariateShapeData<Number2>>
537 &shape_data,
538 Number *values_dofs_in,
539 Number *values,
540 Number *gradients,
541 Number *scratch_data,
542 const unsigned int subface_index,
543 const unsigned int face_direction)
544 {
545 AssertDimension(shape_data.size(), 2);
546
547 const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
548 const int n_rows_n = degree + 1;
549 const int n_rows_t = degree;
550 const ::ndarray<int, 3, 3> dofs_per_direction{
551 {{{n_rows_n, n_rows_t, n_rows_t}},
552 {{n_rows_t, n_rows_n, n_rows_t}},
553 {{n_rows_t, n_rows_t, n_rows_n}}}};
554
555 (void)scratch_data;
556 (void)subface_index;
557 // TODO: This is currently not implemented, but the test
558 // matrix_vector_rt_face_03 apparently works without it -> check
559 // if (subface_index < GeometryInfo<dim - 1>::max_children_per_cell)
560 // DEAL_II_NOT_IMPLEMENTED();
561
562 using Eval = EvaluatorTensorProduct<evaluate_evenodd,
563 dim - 1,
564 (fe_degree > 0 ? fe_degree : 0),
565 n_q_points_1d,
566 Number,
567 Number2>;
568
569 std::array<int, dim> values_dofs_offsets = {};
570 for (unsigned int comp = 0; comp < dim - 1; ++comp)
571 {
572 if (dim == 2)
573 values_dofs_offsets[comp + 1] =
574 values_dofs_offsets[comp] +
575 3 * dofs_per_direction[comp][(face_direction + 1) % dim];
576 else
577 values_dofs_offsets[comp + 1] =
578 values_dofs_offsets[comp] +
579 3 * dofs_per_direction[comp][(face_direction + 1) % dim] *
580 dofs_per_direction[comp][(face_direction + 2) % dim];
581 }
582
583 // Jacobians on faces are reordered to enable simple access with the
584 // regular evaluators; to get the RT Piola transform right, we need to
585 // pass through the values_dofs array in a permuted right order
586 std::array<unsigned int, dim> components;
587 for (unsigned int comp = 0; comp < dim; ++comp)
588 components[comp] = (face_direction + comp + 1) % dim;
589
590 for (const unsigned int comp : components)
591 {
592 Number *values_dofs = values_dofs_in + values_dofs_offsets[comp];
593
594 std::array<int, 2> n_blocks{
595 {dofs_per_direction[comp][(face_direction + 1) % dim],
596 (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] :
597 1)}};
598
599 if constexpr (dim == 3)
600 {
601 EvaluatorTensorProduct<evaluate_evenodd,
602 dim - 1,
603 n_q_points_1d,
604 n_q_points_1d,
605 Number,
606 Number2>
607 eval_g({},
608 shape_data[0].shape_gradients_collocation_eo.data(),
609 {});
610 if (!do_integrate)
611 {
612 EvaluatorTensorProductAnisotropic<dim - 1,
613 fe_degree,
614 n_q_points_1d,
615 true>
616 eval;
617 // Evaluate in 3d
618 if (n_blocks[0] == n_rows_n)
619 {
620 eval.template normal<0>(shape_data[0],
621 values_dofs,
622 values);
623 eval.template tangential<1, 0>(shape_data[1],
624 values,
625 values);
626
627 if (evaluation_flag & EvaluationFlags::gradients)
628 {
629 eval.template normal<0>(shape_data[0],
630 values_dofs +
631 n_blocks[0] * n_blocks[1],
632 scratch_data);
633 eval.template tangential<1, 0, dim>(shape_data[1],
634 scratch_data,
635 gradients + 2);
636 }
637 }
638 else if (n_blocks[1] == n_rows_n)
639 {
640 eval.template normal<1>(shape_data[0],
641 values_dofs,
642 values);
643 eval.template tangential<0, 1>(shape_data[1],
644 values,
645 values);
646
647 if (evaluation_flag & EvaluationFlags::gradients)
648 {
649 eval.template normal<1>(shape_data[0],
650 values_dofs +
651 n_blocks[0] * n_blocks[1],
652 scratch_data);
653 eval.template tangential<0, 1, dim>(shape_data[1],
654 scratch_data,
655 gradients + 2);
656 }
657 }
658 else
659 {
660 Eval eval(shape_data[1].shape_values_eo.data(), {}, {});
661 eval.template values<0, true, false>(values_dofs, values);
662 eval.template values<1, true, false>(values, values);
663 if (evaluation_flag & EvaluationFlags::gradients)
664 {
665 eval.template values<0, true, false>(values_dofs +
666 n_blocks[0] *
667 n_blocks[1],
668 scratch_data);
669 eval.template values<1, true, false, dim>(
670 scratch_data, gradients + 2);
671 }
672 }
673 if (evaluation_flag & EvaluationFlags::gradients)
674 {
675 eval_g.template gradients<0, true, false, dim>(values,
676 gradients);
677 eval_g.template gradients<1, true, false, dim>(values,
678 gradients +
679 1);
680 }
681 }
682 else
683 {
684 EvaluatorTensorProductAnisotropic<dim - 1,
685 fe_degree,
686 n_q_points_1d,
687 false>
688 eval;
689 // Integrate in 3d
690 if (evaluation_flag & EvaluationFlags::gradients)
691 {
692 if (evaluation_flag & EvaluationFlags::values)
693 eval_g.template gradients<0, false, true, dim>(
694 gradients, values);
695 else
696 eval_g.template gradients<0, false, false, dim>(
697 gradients, values);
698 eval_g.template gradients<1, false, true, dim>(gradients +
699 1,
700 values);
701 }
702 if (n_blocks[0] == n_rows_n)
703 {
704 eval.template tangential<1, 0>(shape_data[1],
705 values,
706 values);
707 eval.template normal<0>(shape_data[0],
708 values,
709 values_dofs);
710
711 if (evaluation_flag & EvaluationFlags::gradients)
712 {
713 eval.template tangential<1, 0, dim>(shape_data[1],
714 gradients + 2,
715 scratch_data);
716 eval.template normal<0>(shape_data[0],
717 scratch_data,
718 values_dofs +
719 n_blocks[0] * n_blocks[1]);
720 }
721 }
722 else if (n_blocks[1] == n_rows_n)
723 {
724 eval.template tangential<0, 1>(shape_data[1],
725 values,
726 values);
727 eval.template normal<1>(shape_data[0],
728 values,
729 values_dofs);
730
731 if (evaluation_flag & EvaluationFlags::gradients)
732 {
733 eval.template tangential<0, 1, dim>(shape_data[1],
734 gradients + 2,
735 scratch_data);
736 eval.template normal<1>(shape_data[0],
737 scratch_data,
738 values_dofs +
739 n_blocks[0] * n_blocks[1]);
740 }
741 }
742 else
743 {
744 Eval eval_iso(shape_data[1].shape_values_eo.data(),
745 {},
746 {});
747 eval_iso.template values<1, false, false>(values, values);
748 eval_iso.template values<0, false, false>(values,
749 values_dofs);
750 if (evaluation_flag & EvaluationFlags::gradients)
751 {
752 eval_iso.template values<1, false, false, dim>(
753 gradients + 2, scratch_data);
754 eval_iso.template values<0, false, false>(
755 scratch_data,
756 values_dofs + n_blocks[0] * n_blocks[1]);
757 }
758 }
759 }
760 }
761 else
762 {
763 using EvalN = EvaluatorTensorProduct<evaluate_evenodd,
764 dim - 1,
765 fe_degree + 1,
766 n_q_points_1d,
767 Number,
768 Number2>;
769 if (!do_integrate)
770 {
771 // Evaluate in 2d
772 if (n_blocks[0] == n_rows_n)
773 {
774 EvalN eval(shape_data[0].shape_values_eo,
775 shape_data[0].shape_gradients_eo,
776 {});
777 eval.template values<0, true, false>(values_dofs, values);
778 if (evaluation_flag & EvaluationFlags::gradients)
779 {
780 eval.template gradients<0, true, false, dim>(
781 values_dofs, gradients);
782 eval.template values<0, true, false, dim>(
783 values_dofs + n_rows_n, gradients + 1);
784 }
785 }
786 else
787 {
788 Eval eval(shape_data[1].shape_values_eo,
789 shape_data[1].shape_gradients_eo,
790 {});
791 eval.template values<0, true, false>(values_dofs, values);
792 if (evaluation_flag & EvaluationFlags::gradients)
793 {
794 eval.template gradients<0, true, false, dim>(
795 values_dofs, gradients);
796 eval.template values<0, true, false, dim>(
797 values_dofs + n_rows_t, gradients + 1);
798 }
799 }
800 }
801 else
802 {
803 // Integrate in 2d
804 if (n_blocks[0] == n_rows_n)
805 {
806 EvalN eval(shape_data[0].shape_values_eo,
807 shape_data[0].shape_gradients_eo,
808 {});
809 if (evaluation_flag & EvaluationFlags::values)
810 eval.template values<0, false, false>(values,
811 values_dofs);
812 if (evaluation_flag & EvaluationFlags::gradients)
813 {
814 if (evaluation_flag & EvaluationFlags::values)
815 eval.template gradients<0, false, true, dim>(
816 gradients, values_dofs);
817 else
818 eval.template gradients<0, false, false, dim>(
819 gradients, values_dofs);
820 eval.template values<0, false, false, dim>(
821 gradients + 1, values_dofs + n_rows_n);
822 }
823 }
824 else
825 {
826 Eval eval(shape_data[1].shape_values_eo,
827 shape_data[1].shape_gradients_eo,
828 {});
829 if (evaluation_flag & EvaluationFlags::values)
830 eval.template values<0, false, false>(values,
831 values_dofs);
832 if (evaluation_flag & EvaluationFlags::gradients)
833 {
834 if (evaluation_flag & EvaluationFlags::values)
835 eval.template gradients<0, false, true, dim>(
836 gradients, values_dofs);
837 else
838 eval.template gradients<0, false, false, dim>(
839 gradients, values_dofs);
840 eval.template values<0, false, false, dim>(
841 gradients + 1, values_dofs + n_rows_t);
842 }
843 }
844 }
845 }
846 values += Utilities::pow(n_q_points_1d, dim - 1);
847 gradients += dim * Utilities::pow(n_q_points_1d, dim - 1);
848 }
849 }
850 };
851
852
853
854 template <int dim, int fe_degree, typename Number>
855 struct FEFaceNormalEvaluationImpl
856 {
857 using Number2 =
858 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
859
860 template <bool do_evaluate, bool add_into_output>
861 static void
862 interpolate(const unsigned int n_components,
863 const EvaluationFlags::EvaluationFlags flags,
864 const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
865 const Number *input,
866 Number *output,
867 const unsigned int face_no)
868 {
869 Assert(static_cast<unsigned int>(fe_degree) ==
870 shape_info.data.front().fe_degree ||
871 fe_degree == -1,
872 ExcInternalError());
873 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
874 interpolate_raviart_thomas<do_evaluate, add_into_output>(
875 n_components, input, output, flags, face_no, shape_info);
876 else
877 interpolate_generic<do_evaluate, add_into_output>(
878 n_components,
879 input,
880 output,
881 flags,
882 face_no,
883 shape_info.data.front().fe_degree + 1,
884 shape_info.data.front().shape_data_on_face,
885 shape_info.dofs_per_component_on_cell,
886 3 * shape_info.dofs_per_component_on_face);
887 }
888
892 template <bool do_evaluate, bool add_into_output>
893 static void
894 interpolate_quadrature(
895 const unsigned int n_components,
896 const EvaluationFlags::EvaluationFlags flags,
897 const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
898 const Number *input,
899 Number *output,
900 const unsigned int face_no)
901 {
902 Assert(static_cast<unsigned int>(fe_degree + 1) ==
903 shape_info.data.front().n_q_points_1d ||
904 fe_degree == -1,
905 ExcInternalError());
906
907 interpolate_generic<do_evaluate, add_into_output>(
908 n_components,
909 input,
910 output,
911 flags,
912 face_no,
913 shape_info.data.front().quadrature.size(),
914 shape_info.data.front().quadrature_data_on_face,
915 shape_info.n_q_points,
916 shape_info.n_q_points_face);
917 }
918
919 private:
920 template <bool do_evaluate, bool add_into_output, int face_direction = 0>
921 static void
922 interpolate_generic(const unsigned int n_components,
923 const Number *input,
924 Number *output,
925 const EvaluationFlags::EvaluationFlags flag,
926 const unsigned int face_no,
927 const unsigned int n_points_1d,
928 const std::array<AlignedVector<Number2>, 2> &shape_data,
929 const unsigned int dofs_per_component_on_cell,
930 const unsigned int dofs_per_component_on_face)
931 {
932 if (face_direction == face_no / 2)
933 {
934 constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
935
936 const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
937 const int stride = Utilities::pow(n_rows, face_direction);
938 const std::array<int, 2> n_blocks{
939 {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
940 std::array<int, 2> steps;
941 if constexpr (face_direction == 0)
942 steps = {{n_rows, 0}};
943 else if constexpr (face_direction == 1 && dim == 2)
944 steps = {{1, 0}};
945 else if constexpr (face_direction == 1)
946 // in 3d, the coordinate system is zx, not xz -> switch indices
947 steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
948 else if constexpr (face_direction == 2)
949 steps = {{1, 0}};
950
951 for (unsigned int c = 0; c < n_components; ++c)
952 {
953 if (flag & EvaluationFlags::hessians)
954 interpolate_to_face<fe_degree + 1,
955 stride_,
956 do_evaluate,
957 add_into_output,
958 2>(shape_data[face_no % 2].begin(),
959 n_blocks,
960 steps,
961 input,
962 output,
963 n_rows,
964 stride);
965 else if (flag & EvaluationFlags::gradients)
966 interpolate_to_face<fe_degree + 1,
967 stride_,
968 do_evaluate,
969 add_into_output,
970 1>(shape_data[face_no % 2].begin(),
971 n_blocks,
972 steps,
973 input,
974 output,
975 n_rows,
976 stride);
977 else
978 interpolate_to_face<fe_degree + 1,
979 stride_,
980 do_evaluate,
981 add_into_output,
982 0>(shape_data[face_no % 2].begin(),
983 n_blocks,
984 steps,
985 input,
986 output,
987 n_rows,
988 stride);
989 if (do_evaluate)
990 {
991 input += dofs_per_component_on_cell;
992 output += dofs_per_component_on_face;
993 }
994 else
995 {
996 output += dofs_per_component_on_cell;
997 input += dofs_per_component_on_face;
998 }
999 }
1000 }
1001 else if (face_direction < dim)
1002 {
1003 interpolate_generic<do_evaluate,
1004 add_into_output,
1005 std::min(face_direction + 1, dim - 1)>(
1006 n_components,
1007 input,
1008 output,
1009 flag,
1010 face_no,
1011 n_points_1d,
1012 shape_data,
1013 dofs_per_component_on_cell,
1014 dofs_per_component_on_face);
1015 }
1016 }
1017
1018 template <bool do_evaluate,
1019 bool add_into_output,
1020 int face_direction = 0,
1021 int max_derivative = 0>
1022 static void
1023 interpolate_raviart_thomas(
1024 const unsigned int n_components,
1025 const Number *input,
1026 Number *output,
1027 const EvaluationFlags::EvaluationFlags flag,
1028 const unsigned int face_no,
1029 const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info)
1030 {
1031 if (dim == 1)
1032 {
1033 // This should never happen since the FE_RaviartThomasNodal is not
1034 // defined for dim = 1. It prevents compiler warnings of infinite
1035 // recursion.
1036 DEAL_II_ASSERT_UNREACHABLE();
1037 return;
1038 }
1039
1040 bool increase_max_der = false;
1041 if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
1042 (flag & EvaluationFlags::gradients && max_derivative < 1))
1043 increase_max_der = true;
1044
1045 if (face_direction == face_no / 2 && !increase_max_der)
1046 {
1047 constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
1048 constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
1049 constexpr int stride2 = fe_degree * (fe_degree + 1);
1050
1051 const int degree =
1052 fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
1053 const int n_rows_n = degree + 1;
1054 const int n_rows_t = degree;
1055
1056 std::array<int, 3> strides{{1, 1, 1}};
1057 if (face_direction > 0)
1058 {
1059 strides[0] =
1060 n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
1061 strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
1062 strides[2] = Utilities::pow(n_rows_t, face_direction);
1063 }
1064 const ::ndarray<int, 3, 3> dofs_per_direction{
1065 {{{n_rows_n, n_rows_t, n_rows_t}},
1066 {{n_rows_t, n_rows_n, n_rows_t}},
1067 {{n_rows_t, n_rows_t, n_rows_n}}}};
1068
1069 std::array<int, 2> steps, n_blocks;
1070
1071 if constexpr (face_direction == 0)
1072 steps = {{degree + (face_direction == 0), 0}};
1073 else if constexpr (face_direction == 1 && dim == 2)
1074 steps = {{1, 0}};
1075 else if constexpr (face_direction == 1)
1076 // in 3d, the coordinate system is zx, not xz -> switch indices
1077 steps = {
1078 {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1079 else if constexpr (face_direction == 2)
1080 steps = {{1, 0}};
1081
1082 n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
1083 n_blocks[1] =
1084 dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
1085
1086 interpolate_to_face<
1087 (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
1088 ((face_direction < 2) ? stride1 : stride2),
1089 do_evaluate,
1090 add_into_output,
1091 max_derivative>(shape_info.data[face_direction != 0]
1092 .shape_data_on_face[face_no % 2]
1093 .begin(),
1094 n_blocks,
1095 steps,
1096 input,
1097 output,
1098 degree + (face_direction == 0),
1099 strides[0]);
1100
1101 if (do_evaluate)
1102 {
1103 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1104 output += 3 * n_blocks[0] * n_blocks[1];
1105 }
1106 else
1107 {
1108 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1109 input += 3 * n_blocks[0] * n_blocks[1];
1110 }
1111
1112 // must only change steps only for face direction 0
1113 if constexpr (face_direction == 0)
1114 steps = {{degree, 0}};
1115
1116 n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
1117 n_blocks[1] =
1118 dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
1119
1120 interpolate_to_face<
1121 (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
1122 ((face_direction < 2) ? stride0 : stride2),
1123 do_evaluate,
1124 add_into_output,
1125 max_derivative>(shape_info.data[face_direction != 1]
1126 .shape_data_on_face[face_no % 2]
1127 .begin(),
1128 n_blocks,
1129 steps,
1130 input,
1131 output,
1132 degree + (face_direction == 1),
1133 strides[1]);
1134
1135 if constexpr (dim > 2)
1136 {
1137 if (do_evaluate)
1138 {
1139 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1140 output += 3 * n_blocks[0] * n_blocks[1];
1141 }
1142 else
1143 {
1144 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1145 input += 3 * n_blocks[0] * n_blocks[1];
1146 }
1147
1148 if constexpr (face_direction == 0)
1149 steps = {{degree, 0}};
1150 else if constexpr (face_direction == 1)
1151 // in 3d, the coordinate system is zx, not xz -> switch indices
1152 steps = {
1153 {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1154 else if constexpr (face_direction == 2)
1155 steps = {{1, 0}};
1156
1157 n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
1158 n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
1159
1160 interpolate_to_face<
1161 (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
1162 stride0,
1163 do_evaluate,
1164 add_into_output,
1165 max_derivative>(shape_info.data[face_direction != 2]
1166 .shape_data_on_face[face_no % 2]
1167 .begin(),
1168 n_blocks,
1169 steps,
1170 input,
1171 output,
1172 degree + (face_direction == 2),
1173 strides[2]);
1174 }
1175 }
1176 else if (face_direction == face_no / 2)
1177 {
1178 // Only increase max_derivative
1179 interpolate_raviart_thomas<do_evaluate,
1180 add_into_output,
1181 face_direction,
1182 std::min(max_derivative + 1, 2)>(
1183 n_components, input, output, flag, face_no, shape_info);
1184 }
1185 else if (face_direction < dim)
1186 {
1187 if (increase_max_der)
1188 {
1189 interpolate_raviart_thomas<do_evaluate,
1190 add_into_output,
1191 std::min(face_direction + 1, dim - 1),
1192 std::min(max_derivative + 1, 2)>(
1193 n_components, input, output, flag, face_no, shape_info);
1194 }
1195 else
1196 {
1197 interpolate_raviart_thomas<do_evaluate,
1198 add_into_output,
1199 std::min(face_direction + 1, dim - 1),
1200 max_derivative>(
1201 n_components, input, output, flag, face_no, shape_info);
1202 }
1203 }
1204 }
1205 };
1206
1207
1208
1209 // internal helper function for reading data; base version of different types
1210 template <typename VectorizedArrayType, typename Number2>
1211 void
1212 do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
1213 {
1214 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1215 dst[v] = src_ptr[v];
1216 }
1217
1218
1219
1220 // internal helper function for reading data; specialized version where we
1221 // can use a dedicated load function
1222 template <typename Number, std::size_t width>
1223 void
1224 do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
1225 {
1226 dst.load(src_ptr);
1227 }
1228
1229
1230
1231 // internal helper function for reading data; base version of different types
1232 template <typename VectorizedArrayType, typename Number2>
1233 void
1234 do_vectorized_gather(const Number2 *src_ptr,
1235 const unsigned int *indices,
1236 VectorizedArrayType &dst)
1237 {
1238 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1239 dst[v] = src_ptr[indices[v]];
1240 }
1241
1242
1243
1244 // internal helper function for reading data; specialized version where we
1245 // can use a dedicated gather function
1246 template <typename Number, std::size_t width>
1247 void
1248 do_vectorized_gather(const Number *src_ptr,
1249 const unsigned int *indices,
1250 VectorizedArray<Number, width> &dst)
1251 {
1252 dst.gather(src_ptr, indices);
1253 }
1254
1255
1256
1257 // internal helper function for reading data; base version of different types
1258 template <typename VectorizedArrayType, typename Number2>
1259 void
1260 do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
1261 {
1262 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1263 dst_ptr[v] += src[v];
1264 }
1265
1266
1267
1268 // internal helper function for reading data; specialized version where we
1269 // can use a dedicated load function
1270 template <typename Number, std::size_t width>
1271 void
1272 do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
1273 {
1274 VectorizedArray<Number, width> tmp;
1275 tmp.load(dst_ptr);
1276 (tmp + src).store(dst_ptr);
1277 }
1278
1279
1280
1281 // internal helper function for reading data; base version of different types
1282 template <typename VectorizedArrayType, typename Number2>
1283 void
1284 do_vectorized_scatter_add(const VectorizedArrayType src,
1285 const unsigned int *indices,
1286 Number2 *dst_ptr)
1287 {
1288 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1289 dst_ptr[indices[v]] += src[v];
1290 }
1291
1292
1293
1294 // internal helper function for reading data; specialized version where we
1295 // can use a dedicated gather function
1296 template <typename Number, std::size_t width>
1297 void
1298 do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
1299 const unsigned int *indices,
1300 Number *dst_ptr)
1301 {
1302#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
1303 for (unsigned int v = 0; v < width; ++v)
1304 dst_ptr[indices[v]] += src[v];
1305#else
1306 VectorizedArray<Number, width> tmp;
1307 tmp.gather(dst_ptr, indices);
1308 (tmp + src).scatter(indices, dst_ptr);
1309#endif
1310 }
1311
1312
1313
1314 template <typename Number>
1315 void
1316 adjust_for_face_orientation(const unsigned int dim,
1317 const unsigned int n_components,
1318 const EvaluationFlags::EvaluationFlags flag,
1319 const unsigned int *orientation,
1320 const bool integrate,
1321 const std::size_t n_q_points,
1322 Number *tmp_values,
1323 Number *values_quad,
1324 Number *gradients_quad,
1325 Number *hessians_quad)
1326 {
1327 for (unsigned int c = 0; c < n_components; ++c)
1328 {
1329 if (flag & EvaluationFlags::values)
1330 {
1331 if (integrate)
1332 for (unsigned int q = 0; q < n_q_points; ++q)
1333 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
1334 else
1335 for (unsigned int q = 0; q < n_q_points; ++q)
1336 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
1337 for (unsigned int q = 0; q < n_q_points; ++q)
1338 values_quad[c * n_q_points + q] = tmp_values[q];
1339 }
1340 if (flag & EvaluationFlags::gradients)
1341 for (unsigned int d = 0; d < dim; ++d)
1342 {
1343 if (integrate)
1344 for (unsigned int q = 0; q < n_q_points; ++q)
1345 tmp_values[q] =
1346 gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
1347 else
1348 for (unsigned int q = 0; q < n_q_points; ++q)
1349 tmp_values[orientation[q]] =
1350 gradients_quad[(c * n_q_points + q) * dim + d];
1351 for (unsigned int q = 0; q < n_q_points; ++q)
1352 gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
1353 }
1354 if (flag & EvaluationFlags::hessians)
1355 {
1356 const unsigned int hdim = (dim * (dim + 1)) / 2;
1357 for (unsigned int d = 0; d < hdim; ++d)
1358 {
1359 if (integrate)
1360 for (unsigned int q = 0; q < n_q_points; ++q)
1361 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1362 orientation[q]];
1363 else
1364 for (unsigned int q = 0; q < n_q_points; ++q)
1365 tmp_values[orientation[q]] =
1366 hessians_quad[(c * hdim + d) * n_q_points + q];
1367 for (unsigned int q = 0; q < n_q_points; ++q)
1368 hessians_quad[(c * hdim + d) * n_q_points + q] =
1369 tmp_values[q];
1370 }
1371 }
1372 }
1373 }
1374
1375
1376
1377 template <typename Number, typename VectorizedArrayType>
1378 void
1379 adjust_for_face_orientation_per_lane(
1380 const unsigned int dim,
1381 const unsigned int n_components,
1382 const unsigned int v,
1383 const EvaluationFlags::EvaluationFlags flag,
1384 const unsigned int *orientation,
1385 const bool integrate,
1386 const std::size_t n_q_points,
1387 Number *tmp_values,
1388 VectorizedArrayType *values_quad,
1389 VectorizedArrayType *gradients_quad = nullptr,
1390 VectorizedArrayType *hessians_quad = nullptr)
1391 {
1392 for (unsigned int c = 0; c < n_components; ++c)
1393 {
1394 if (flag & EvaluationFlags::values)
1395 {
1396 if (integrate)
1397 for (unsigned int q = 0; q < n_q_points; ++q)
1398 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
1399 else
1400 for (unsigned int q = 0; q < n_q_points; ++q)
1401 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
1402 for (unsigned int q = 0; q < n_q_points; ++q)
1403 values_quad[c * n_q_points + q][v] = tmp_values[q];
1404 }
1405 if (flag & EvaluationFlags::gradients)
1406 for (unsigned int d = 0; d < dim; ++d)
1407 {
1408 Assert(gradients_quad != nullptr, ExcInternalError());
1409 if (integrate)
1410 for (unsigned int q = 0; q < n_q_points; ++q)
1411 tmp_values[q] =
1412 gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
1413 [v];
1414 else
1415 for (unsigned int q = 0; q < n_q_points; ++q)
1416 tmp_values[orientation[q]] =
1417 gradients_quad[(c * n_q_points + q) * dim + d][v];
1418 for (unsigned int q = 0; q < n_q_points; ++q)
1419 gradients_quad[(c * n_q_points + q) * dim + d][v] =
1420 tmp_values[q];
1421 }
1422 if (flag & EvaluationFlags::hessians)
1423 {
1424 Assert(hessians_quad != nullptr, ExcInternalError());
1425 const unsigned int hdim = (dim * (dim + 1)) / 2;
1426 for (unsigned int d = 0; d < hdim; ++d)
1427 {
1428 if (integrate)
1429 for (unsigned int q = 0; q < n_q_points; ++q)
1430 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1431 orientation[q]][v];
1432 else
1433 for (unsigned int q = 0; q < n_q_points; ++q)
1434 tmp_values[orientation[q]] =
1435 hessians_quad[(c * hdim + d) * n_q_points + q][v];
1436 for (unsigned int q = 0; q < n_q_points; ++q)
1437 hessians_quad[(c * hdim + d) * n_q_points + q][v] =
1438 tmp_values[q];
1439 }
1440 }
1441 }
1442 }
1443
1444
1445
1446 template <int dim, typename Number>
1447 struct FEFaceEvaluationImplEvaluateSelector
1448 {
1449 static bool
1450 evaluate_tensor_none(const unsigned int n_components,
1451 const EvaluationFlags::EvaluationFlags evaluation_flag,
1452 const Number *values_dofs,
1453 FEEvaluationData<dim, Number, true> &fe_eval)
1454 {
1455 const auto &shape_info = fe_eval.get_shape_info();
1456 const auto &shape_data = shape_info.data.front();
1457 using Number2 =
1458 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
1459
1460 Assert((fe_eval.get_dof_access_index() ==
1461 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1462 fe_eval.is_interior_face() == false) == false,
1463 ExcNotImplemented());
1464
1465 const unsigned int face_no = fe_eval.get_face_no();
1466 const unsigned int face_orientation = fe_eval.get_face_orientation();
1467 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1468 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1469
1470 if (evaluation_flag & EvaluationFlags::values)
1471 {
1472 const auto *const shape_values =
1473 &shape_data.shape_values_face(face_no, face_orientation, 0);
1474
1475 auto *out = fe_eval.begin_values();
1476 auto *in = values_dofs;
1477
1478 for (unsigned int c = 0; c < n_components; c += 3)
1479 {
1480 if (c + 1 == n_components)
1481 apply_matrix_vector_product<evaluate_general,
1482 EvaluatorQuantity::value,
1483 /*transpose_matrix*/ true,
1484 /*add*/ false,
1485 /*consider_strides*/ false,
1486 Number,
1487 Number2,
1488 /*n_components*/ 1>(
1489 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1490 else if (c + 2 == n_components)
1491 apply_matrix_vector_product<evaluate_general,
1492 EvaluatorQuantity::value,
1493 /*transpose_matrix*/ true,
1494 /*add*/ false,
1495 /*consider_strides*/ false,
1496 Number,
1497 Number2,
1498 /*n_components*/ 2>(
1499 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1500 else
1501 apply_matrix_vector_product<evaluate_general,
1502 EvaluatorQuantity::value,
1503 /*transpose_matrix*/ true,
1504 /*add*/ false,
1505 /*consider_strides*/ false,
1506 Number,
1507 Number2,
1508 /*n_components*/ 3>(
1509 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1510
1511 out += 3 * n_q_points;
1512 in += 3 * n_dofs;
1513 }
1514 }
1515
1516 if (evaluation_flag & EvaluationFlags::gradients)
1517 {
1518 auto *out = fe_eval.begin_gradients();
1519 const auto *in = values_dofs;
1520
1521 const auto *const shape_gradients =
1522 &shape_data.shape_gradients_face(face_no, face_orientation, 0);
1523
1524 for (unsigned int c = 0; c < n_components; c += 3)
1525 {
1526 if (c + 1 == n_components)
1527 apply_matrix_vector_product<evaluate_general,
1528 EvaluatorQuantity::value,
1529 /*transpose_matrix*/ true,
1530 /*add*/ false,
1531 /*consider_strides*/ false,
1532 Number,
1533 Number2,
1534 /*n_components*/ 1>(
1535 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1536 else if (c + 2 == n_components)
1537 apply_matrix_vector_product<evaluate_general,
1538 EvaluatorQuantity::value,
1539 /*transpose_matrix*/ true,
1540 /*add*/ false,
1541 /*consider_strides*/ false,
1542 Number,
1543 Number2,
1544 /*n_components*/ 2>(
1545 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1546 else
1547 apply_matrix_vector_product<evaluate_general,
1548 EvaluatorQuantity::value,
1549 /*transpose_matrix*/ true,
1550 /*add*/ false,
1551 /*consider_strides*/ false,
1552 Number,
1553 Number2,
1554 /*n_components*/ 3>(
1555 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1556 out += 3 * n_q_points * dim;
1557 in += 3 * n_dofs;
1558 }
1559 }
1560
1561 Assert(!(evaluation_flag & EvaluationFlags::hessians),
1562 ExcNotImplemented());
1563
1564 return true;
1565 }
1566
1567 template <int fe_degree>
1568#ifndef DEBUG
1569 DEAL_II_ALWAYS_INLINE
1570#endif
1571 static void
1572 project_to_face(const unsigned int n_components,
1573 const EvaluationFlags::EvaluationFlags evaluation_flag,
1574 const Number *values_dofs,
1575 FEEvaluationData<dim, Number, true> &fe_eval,
1576 const bool use_vectorization,
1577 Number *temp,
1578 Number *scratch_data)
1579 {
1580 const auto &shape_info = fe_eval.get_shape_info();
1581
1582 if (use_vectorization == false)
1583 {
1584 const auto &shape_data = shape_info.data.front();
1585
1586 const unsigned int dofs_per_comp_face =
1587 fe_degree > -1 ?
1588 Utilities::pow(fe_degree + 1, dim - 1) :
1589 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1590 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
1591
1592 for (unsigned int v = 0; v < Number::size(); ++v)
1593 {
1594 // the loop breaks once an invalid_unsigned_int is hit for
1595 // all cases except the exterior faces in the ECL loop (where
1596 // some faces might be at the boundaries but others not)
1597 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1598 {
1599 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1600 temp[i][v] = 0;
1601 continue;
1602 }
1603
1604 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
1605 template interpolate<true, false>(n_components,
1606 evaluation_flag,
1607 shape_info,
1608 values_dofs,
1609 scratch_data,
1610 fe_eval.get_face_no(v));
1611
1612 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1613 temp[i][v] = scratch_data[i][v];
1614 }
1615 }
1616 else
1617 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
1618 template interpolate<true, false>(n_components,
1619 evaluation_flag,
1620 shape_info,
1621 values_dofs,
1622 temp,
1623 fe_eval.get_face_no());
1624 }
1625
1626
1627 template <int fe_degree, int n_q_points_1d>
1628#ifndef DEBUG
1629 DEAL_II_ALWAYS_INLINE
1630#endif
1631 static void
1632 evaluate_in_face(const unsigned int n_components,
1633 const EvaluationFlags::EvaluationFlags evaluation_flag,
1634 FEEvaluationData<dim, Number, true> &fe_eval,
1635 Number *temp,
1636 Number *scratch_data)
1637 {
1638 const auto &shape_info = fe_eval.get_shape_info();
1639 const auto &shape_data = shape_info.data.front();
1640
1641 const unsigned int subface_index = fe_eval.get_subface_index();
1642 constexpr unsigned int n_q_points_1d_actual =
1643 fe_degree > -1 ? n_q_points_1d : 0;
1644
1645 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
1646 {
1647 FEFaceEvaluationImplRaviartThomas<dim,
1648 fe_degree,
1649 n_q_points_1d_actual,
1650 Number>::
1651 template evaluate_or_integrate_in_face<false>(
1652 evaluation_flag,
1653 fe_eval.get_shape_info().data,
1654 temp,
1655 fe_eval.begin_values(),
1656 fe_eval.begin_gradients(),
1657 scratch_data,
1658 subface_index,
1659 fe_eval.get_face_no() / 2);
1660 }
1661 else if (fe_degree > -1 &&
1662 subface_index >= GeometryInfo<dim>::max_children_per_cell &&
1663 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
1664 FEFaceEvaluationImpl<true,
1665 dim,
1666 fe_degree,
1667 n_q_points_1d_actual,
1668 Number>::evaluate_in_face(n_components,
1669 evaluation_flag,
1670 shape_data,
1671 temp,
1672 fe_eval.begin_values(),
1673 fe_eval
1674 .begin_gradients(),
1675 fe_eval.begin_hessians(),
1676 scratch_data,
1677 subface_index);
1678 else
1679 FEFaceEvaluationImpl<false,
1680 dim,
1681 fe_degree,
1682 n_q_points_1d_actual,
1683 Number>::evaluate_in_face(n_components,
1684 evaluation_flag,
1685 shape_data,
1686 temp,
1687 fe_eval.begin_values(),
1688 fe_eval
1689 .begin_gradients(),
1690 fe_eval.begin_hessians(),
1691 scratch_data,
1692 subface_index);
1693 }
1694
1695#ifndef DEBUG
1696 DEAL_II_ALWAYS_INLINE
1697#endif
1698 static void
1699 adjust_quadrature_for_face_orientation(
1700 const unsigned int n_components,
1701 const EvaluationFlags::EvaluationFlags evaluation_flag,
1702 FEEvaluationData<dim, Number, true> &fe_eval,
1703 const bool use_vectorization,
1704 Number *temp)
1705 {
1706 const auto &shape_info = fe_eval.get_shape_info();
1707
1708 if (use_vectorization == false)
1709 {
1710 for (unsigned int v = 0; v < Number::size(); ++v)
1711 {
1712 // the loop breaks once an invalid_unsigned_int is hit for
1713 // all cases except the exterior faces in the ECL loop (where
1714 // some faces might be at the boundaries but others not)
1715 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1716 continue;
1717
1718 if (fe_eval.get_face_orientation(v) != 0)
1719 adjust_for_face_orientation_per_lane(
1720 dim,
1721 n_components,
1722 v,
1723 evaluation_flag,
1724 &shape_info.face_orientations_quad(
1725 fe_eval.get_face_orientation(v), 0),
1726 false,
1727 shape_info.n_q_points_face,
1728 &temp[0][0],
1729 fe_eval.begin_values(),
1730 fe_eval.begin_gradients(),
1731 fe_eval.begin_hessians());
1732 }
1733 }
1734 else if (fe_eval.get_face_orientation() != 0)
1735 adjust_for_face_orientation(
1736 dim,
1737 n_components,
1738 evaluation_flag,
1739 &shape_info.face_orientations_quad(fe_eval.get_face_orientation(), 0),
1740 false,
1741 shape_info.n_q_points_face,
1742 temp,
1743 fe_eval.begin_values(),
1744 fe_eval.begin_gradients(),
1745 fe_eval.begin_hessians());
1746 }
1747
1748
1749
1750 template <int fe_degree, int n_q_points_1d>
1751 static bool
1752 evaluate_tensor(const unsigned int n_components,
1753 const EvaluationFlags::EvaluationFlags evaluation_flag,
1754 const Number *values_dofs,
1755 FEEvaluationData<dim, Number, true> &fe_eval)
1756 {
1757 const auto &shape_info = fe_eval.get_shape_info();
1758 const auto &shape_data = shape_info.data.front();
1759
1760 const unsigned int dofs_per_comp_face =
1761 fe_degree > -1 ?
1762 Utilities::pow(fe_degree + 1, dim - 1) :
1763 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1764
1765 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1766 // array, so reserve space for all three here
1767 Number *temp = fe_eval.get_scratch_data().begin();
1768 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1769
1770 bool use_vectorization = true;
1771 if (fe_eval.get_dof_access_index() ==
1772 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1773 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1774 for (unsigned int v = 0; v < Number::size(); ++v)
1775 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1776 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1777 use_vectorization = false;
1778
1779 project_to_face<fe_degree>(n_components,
1780 evaluation_flag,
1781 values_dofs,
1782 fe_eval,
1783 use_vectorization,
1784 temp,
1785 scratch_data);
1786
1787 evaluate_in_face<fe_degree, n_q_points_1d>(
1788 n_components, evaluation_flag, fe_eval, temp, scratch_data);
1789
1790 if (dim == 3)
1791 adjust_quadrature_for_face_orientation(
1792 n_components, evaluation_flag, fe_eval, use_vectorization, temp);
1793
1794 return false;
1795 }
1796
1797 template <int fe_degree, int n_q_points_1d>
1798 static bool
1799 run(const unsigned int n_components,
1800 const EvaluationFlags::EvaluationFlags evaluation_flag,
1801 const Number *values_dofs,
1802 FEEvaluationData<dim, Number, true> &fe_eval)
1803 {
1804 const auto &shape_info = fe_eval.get_shape_info();
1805
1806 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
1807 return evaluate_tensor_none(n_components,
1808 evaluation_flag,
1809 values_dofs,
1810 fe_eval);
1811 else
1812 return evaluate_tensor<fe_degree, n_q_points_1d>(n_components,
1813 evaluation_flag,
1814 values_dofs,
1815 fe_eval);
1816 }
1817 };
1818
1819
1820
1821 template <int dim, typename Number>
1822 struct FEFaceEvaluationImplProjectToFaceSelector
1823 {
1824 template <int fe_degree>
1825 static bool
1826 run(const unsigned int n_components,
1827 const EvaluationFlags::EvaluationFlags evaluation_flag,
1828 const Number *values_dofs,
1829 FEEvaluationData<dim, Number, true> &fe_eval)
1830 {
1831 const auto &shape_info = fe_eval.get_shape_info();
1832 const auto &shape_data = shape_info.data.front();
1833
1834 const unsigned int dofs_per_comp_face =
1835 fe_degree > -1 ?
1836 Utilities::pow(fe_degree + 1, dim - 1) :
1837 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1838
1839 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1840 // array, so reserve space for all three here
1841 Number *temp = fe_eval.get_scratch_data().begin();
1842 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1843
1844 bool use_vectorization = true;
1845 if (fe_eval.get_dof_access_index() ==
1846 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1847 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1848 for (unsigned int v = 0; v < Number::size(); ++v)
1849 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1850 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1851 use_vectorization = false;
1852
1853 FEFaceEvaluationImplEvaluateSelector<dim, Number>::
1854 template project_to_face<fe_degree>(n_components,
1855 evaluation_flag,
1856 values_dofs,
1857 fe_eval,
1858 use_vectorization,
1859 temp,
1860 scratch_data);
1861
1862 return false;
1863 }
1864 };
1865
1866
1867
1868 template <int dim, typename Number>
1869 struct FEFaceEvaluationImplEvaluateInFaceSelector
1870 {
1871 template <int fe_degree, int n_q_points_1d>
1872 static bool
1873 run(const unsigned int n_components,
1874 const EvaluationFlags::EvaluationFlags evaluation_flag,
1875 FEEvaluationData<dim, Number, true> &fe_eval)
1876 {
1877 const auto &shape_info = fe_eval.get_shape_info();
1878 const auto &shape_data = shape_info.data.front();
1879
1880 const unsigned int dofs_per_comp_face =
1881 fe_degree > -1 ?
1882 Utilities::pow(fe_degree + 1, dim - 1) :
1883 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1884
1885 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1886 // array, so reserve space for all three here
1887 Number *temp = fe_eval.get_scratch_data().begin();
1888 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1889
1890 FEFaceEvaluationImplEvaluateSelector<dim, Number>::
1891 template evaluate_in_face<fe_degree, n_q_points_1d>(
1892 n_components, evaluation_flag, fe_eval, temp, scratch_data);
1893
1894 return false;
1895 }
1896 };
1897
1898
1899
1900 template <int dim, typename Number>
1901 struct FEFaceEvaluationImplIntegrateSelector
1902 {
1903 static bool
1904 integrate_tensor_none(
1905 const unsigned int n_components,
1906 const EvaluationFlags::EvaluationFlags integration_flag,
1907 Number *values_dofs,
1908 FEEvaluationData<dim, Number, true> &fe_eval,
1909 const bool sum_into_values)
1910 {
1911 const auto &shape_info = fe_eval.get_shape_info();
1912 const auto &shape_data = shape_info.data.front();
1913 using Number2 =
1914 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
1915
1916 Assert((fe_eval.get_dof_access_index() ==
1917 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1918 fe_eval.is_interior_face() == false) == false,
1919 ExcNotImplemented());
1920
1921 const unsigned int face_no = fe_eval.get_face_no();
1922 const unsigned int face_orientation = fe_eval.get_face_orientation();
1923 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1924 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1925
1926
1927 if (integration_flag & EvaluationFlags::values)
1928 {
1929 const auto *const shape_values =
1930 &shape_data.shape_values_face(face_no, face_orientation, 0);
1931
1932 auto *in = fe_eval.begin_values();
1933 auto *out = values_dofs;
1934
1935 for (unsigned int c = 0; c < n_components; c += 3)
1936 {
1937 if (sum_into_values)
1938 {
1939 if (c + 1 == n_components)
1940 apply_matrix_vector_product<evaluate_general,
1941 EvaluatorQuantity::value,
1942 /*transpose_matrix*/ false,
1943 /*add*/ true,
1944 /*consider_strides*/ false,
1945 Number,
1946 Number2,
1947 /*n_components*/ 1>(
1948 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1949 else if (c + 2 == n_components)
1950 apply_matrix_vector_product<evaluate_general,
1951 EvaluatorQuantity::value,
1952 /*transpose_matrix*/ false,
1953 /*add*/ true,
1954 /*consider_strides*/ false,
1955 Number,
1956 Number2,
1957 /*n_components*/ 2>(
1958 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1959 else
1960 apply_matrix_vector_product<evaluate_general,
1961 EvaluatorQuantity::value,
1962 /*transpose_matrix*/ false,
1963 /*add*/ true,
1964 /*consider_strides*/ false,
1965 Number,
1966 Number2,
1967 /*n_components*/ 3>(
1968 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1969 }
1970 else
1971 {
1972 if (c + 1 == n_components)
1973 apply_matrix_vector_product<evaluate_general,
1974 EvaluatorQuantity::value,
1975 /*transpose_matrix*/ false,
1976 /*add*/ false,
1977 /*consider_strides*/ false,
1978 Number,
1979 Number2,
1980 /*n_components*/ 1>(
1981 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1982 else if (c + 2 == n_components)
1983 apply_matrix_vector_product<evaluate_general,
1984 EvaluatorQuantity::value,
1985 /*transpose_matrix*/ false,
1986 /*add*/ false,
1987 /*consider_strides*/ false,
1988 Number,
1989 Number2,
1990 /*n_components*/ 2>(
1991 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1992 else
1993 apply_matrix_vector_product<evaluate_general,
1994 EvaluatorQuantity::value,
1995 /*transpose_matrix*/ false,
1996 /*add*/ false,
1997 /*consider_strides*/ false,
1998 Number,
1999 Number2,
2000 /*n_components*/ 3>(
2001 shape_values, in, out, n_dofs, n_q_points, 1, 1);
2002 }
2003 in += 3 * n_q_points;
2004 out += 3 * n_dofs;
2005 }
2006 }
2007
2008 if (integration_flag & EvaluationFlags::gradients)
2009 {
2010 auto *in = fe_eval.begin_gradients();
2011 auto *out = values_dofs;
2012
2013 const auto *const shape_gradients =
2014 &shape_data.shape_gradients_face(face_no, face_orientation, 0);
2015
2016 for (unsigned int c = 0; c < n_components; ++c)
2017 {
2018 if (!sum_into_values &&
2019 !(integration_flag & EvaluationFlags::values))
2020 {
2021 if (c + 1 == n_components)
2022 apply_matrix_vector_product<evaluate_general,
2023 EvaluatorQuantity::value,
2024 /*transpose_matrix*/ false,
2025 /*add*/ false,
2026 /*consider_strides*/ false,
2027 Number,
2028 Number2,
2029 /*n_components*/ 1>(
2030 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2031 else if (c + 2 == n_components)
2032 apply_matrix_vector_product<evaluate_general,
2033 EvaluatorQuantity::value,
2034 /*transpose_matrix*/ false,
2035 /*add*/ false,
2036 /*consider_strides*/ false,
2037 Number,
2038 Number2,
2039 /*n_components*/ 2>(
2040 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2041 else
2042 apply_matrix_vector_product<evaluate_general,
2043 EvaluatorQuantity::value,
2044 /*transpose_matrix*/ false,
2045 /*add*/ false,
2046 /*consider_strides*/ false,
2047 Number,
2048 Number2,
2049 /*n_components*/ 3>(
2050 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2051 }
2052 else
2053 {
2054 if (c + 1 == n_components)
2055 apply_matrix_vector_product<evaluate_general,
2056 EvaluatorQuantity::value,
2057 /*transpose_matrix*/ false,
2058 /*add*/ true,
2059 /*consider_strides*/ false,
2060 Number,
2061 Number2,
2062 /*n_components*/ 1>(
2063 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2064 else if (c + 2 == n_components)
2065 apply_matrix_vector_product<evaluate_general,
2066 EvaluatorQuantity::value,
2067 /*transpose_matrix*/ false,
2068 /*add*/ true,
2069 /*consider_strides*/ false,
2070 Number,
2071 Number2,
2072 /*n_components*/ 2>(
2073 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2074 else
2075 apply_matrix_vector_product<evaluate_general,
2076 EvaluatorQuantity::value,
2077 /*transpose_matrix*/ false,
2078 /*add*/ true,
2079 /*consider_strides*/ false,
2080 Number,
2081 Number2,
2082 /*n_components*/ 3>(
2083 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2084 }
2085 in += 3 * n_q_points * dim;
2086 out += 3 * n_dofs;
2087 }
2088 }
2089
2090 Assert(!(integration_flag & EvaluationFlags::hessians),
2091 ExcNotImplemented());
2092
2093 return true;
2094 }
2095
2096#ifndef DEBUG
2097 DEAL_II_ALWAYS_INLINE
2098#endif
2099 static void
2100 adjust_quadrature_for_face_orientation(
2101 const unsigned int n_components,
2102 const EvaluationFlags::EvaluationFlags integration_flag,
2103 FEEvaluationData<dim, Number, true> &fe_eval,
2104 const bool use_vectorization,
2105 Number *temp)
2106 {
2107 const auto &shape_info = fe_eval.get_shape_info();
2108
2109 if (use_vectorization == false)
2110 {
2111 for (unsigned int v = 0; v < Number::size(); ++v)
2112 {
2113 // the loop breaks once an invalid_unsigned_int is hit for
2114 // all cases except the exterior faces in the ECL loop (where
2115 // some faces might be at the boundaries but others not)
2116 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2117 continue;
2118
2119 if (fe_eval.get_face_orientation(v) != 0)
2120 adjust_for_face_orientation_per_lane(
2121 dim,
2122 n_components,
2123 v,
2124 integration_flag,
2125 &fe_eval.get_shape_info().face_orientations_quad(
2126 fe_eval.get_face_orientation(v), 0),
2127 true,
2128 shape_info.n_q_points_face,
2129 &temp[0][0],
2130 fe_eval.begin_values(),
2131 fe_eval.begin_gradients(),
2132 fe_eval.begin_hessians());
2133 }
2134 }
2135 else if (fe_eval.get_face_orientation() != 0)
2136 adjust_for_face_orientation(
2137 dim,
2138 n_components,
2139 integration_flag,
2140 &fe_eval.get_shape_info().face_orientations_quad(
2141 fe_eval.get_face_orientation(), 0),
2142 true,
2143 shape_info.n_q_points_face,
2144 temp,
2145 fe_eval.begin_values(),
2146 fe_eval.begin_gradients(),
2147 fe_eval.begin_hessians());
2148 }
2149
2150 template <int fe_degree, int n_q_points_1d>
2151#ifndef DEBUG
2152 DEAL_II_ALWAYS_INLINE
2153#endif
2154 static void
2155 integrate_in_face(const unsigned int n_components,
2156 const EvaluationFlags::EvaluationFlags integration_flag,
2157 FEEvaluationData<dim, Number, true> &fe_eval,
2158 Number *temp,
2159 Number *scratch_data)
2160 {
2161 const auto &shape_info = fe_eval.get_shape_info();
2162 const auto &shape_data = shape_info.data.front();
2163
2164 const unsigned int n_q_points_1d_actual =
2165 fe_degree > -1 ? n_q_points_1d : 0;
2166 const unsigned int subface_index = fe_eval.get_subface_index();
2167
2168 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
2169 {
2170 FEFaceEvaluationImplRaviartThomas<dim,
2171 fe_degree,
2172 n_q_points_1d_actual,
2173 Number>::
2174 template evaluate_or_integrate_in_face<true>(
2175 integration_flag,
2176 fe_eval.get_shape_info().data,
2177 temp,
2178 fe_eval.begin_values(),
2179 fe_eval.begin_gradients(),
2180 scratch_data,
2181 subface_index,
2182 fe_eval.get_face_no() / 2);
2183 }
2184 else if (fe_degree > -1 &&
2185 fe_eval.get_subface_index() >=
2186 GeometryInfo<dim - 1>::max_children_per_cell &&
2187 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
2188 FEFaceEvaluationImpl<
2189 true,
2190 dim,
2191 fe_degree,
2192 n_q_points_1d_actual,
2193 Number>::integrate_in_face(n_components,
2194 integration_flag,
2195 shape_data,
2196 temp,
2197 fe_eval.begin_values(),
2198 fe_eval.begin_gradients(),
2199 fe_eval.begin_hessians(),
2200 scratch_data,
2201 subface_index);
2202 else
2203 FEFaceEvaluationImpl<
2204 false,
2205 dim,
2206 fe_degree,
2207 n_q_points_1d_actual,
2208 Number>::integrate_in_face(n_components,
2209 integration_flag,
2210 shape_data,
2211 temp,
2212 fe_eval.begin_values(),
2213 fe_eval.begin_gradients(),
2214 fe_eval.begin_hessians(),
2215 scratch_data,
2216 subface_index);
2217 }
2218
2219 template <int fe_degree>
2220#ifndef DEBUG
2221 DEAL_II_ALWAYS_INLINE
2222#endif
2223 static void
2224 collect_from_face(const unsigned int n_components,
2225 const EvaluationFlags::EvaluationFlags integration_flag,
2226 Number *values_dofs,
2227 FEEvaluationData<dim, Number, true> &fe_eval,
2228 const bool use_vectorization,
2229 const Number *temp,
2230 Number *scratch_data,
2231 const bool sum_into_values)
2232 {
2233 const auto &shape_info = fe_eval.get_shape_info();
2234 const auto &shape_data = shape_info.data.front();
2235
2236 const unsigned int dofs_per_comp_face =
2237 fe_degree > -1 ?
2238 Utilities::pow(fe_degree + 1, dim - 1) :
2239 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2240 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
2241
2242 if (use_vectorization == false)
2243 {
2244 for (unsigned int v = 0; v < Number::size(); ++v)
2245 {
2246 // the loop breaks once an invalid_unsigned_int is hit for
2247 // all cases except the exterior faces in the ECL loop (where
2248 // some faces might be at the boundaries but others not)
2249 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2250 continue;
2251
2252 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
2253 template interpolate<false, false>(n_components,
2254 integration_flag,
2255 shape_info,
2256 temp,
2257 scratch_data,
2258 fe_eval.get_face_no(v));
2259
2260 if (sum_into_values)
2261 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2262 values_dofs[i][v] += scratch_data[i][v];
2263 else
2264 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2265 values_dofs[i][v] = scratch_data[i][v];
2266 }
2267 }
2268 else
2269 {
2270 if (sum_into_values)
2271 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
2272 template interpolate<false, true>(n_components,
2273 integration_flag,
2274 shape_info,
2275 temp,
2276 values_dofs,
2277 fe_eval.get_face_no());
2278 else
2279 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
2280 template interpolate<false, false>(n_components,
2281 integration_flag,
2282 shape_info,
2283 temp,
2284 values_dofs,
2285 fe_eval.get_face_no());
2286 }
2287 }
2288
2289 template <int fe_degree, int n_q_points_1d>
2290 static bool
2291 integrate_tensor(const unsigned int n_components,
2292 const EvaluationFlags::EvaluationFlags integration_flag,
2293 Number *values_dofs,
2294 FEEvaluationData<dim, Number, true> &fe_eval,
2295 const bool sum_into_values)
2296 {
2297 const auto &shape_info = fe_eval.get_shape_info();
2298 const auto &shape_data = shape_info.data.front();
2299
2300 const unsigned int dofs_per_comp_face =
2301 fe_degree > -1 ?
2302 Utilities::pow(fe_degree + 1, dim - 1) :
2303 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2304
2305 Number *temp = fe_eval.get_scratch_data().begin();
2306 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2307
2308 bool use_vectorization = true;
2309
2310 if (fe_eval.get_dof_access_index() ==
2311 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
2312 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2313 use_vectorization =
2314 fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
2315 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2316 fe_eval.get_cell_ids().end(),
2317 [&](const auto &v) {
2318 return v == fe_eval.get_cell_ids()[0] ||
2319 v == numbers::invalid_unsigned_int;
2320 });
2321
2322 if (dim == 3)
2323 adjust_quadrature_for_face_orientation(
2324 n_components, integration_flag, fe_eval, use_vectorization, temp);
2325
2326 integrate_in_face<fe_degree, n_q_points_1d>(
2327 n_components, integration_flag, fe_eval, temp, scratch_data);
2328
2329 collect_from_face<fe_degree>(n_components,
2330 integration_flag,
2331 values_dofs,
2332 fe_eval,
2333 use_vectorization,
2334 temp,
2335 scratch_data,
2336 sum_into_values);
2337
2338 return false;
2339 }
2340
2341 template <int fe_degree, int n_q_points_1d>
2342 static bool
2343 run(const unsigned int n_components,
2344 const EvaluationFlags::EvaluationFlags integration_flag,
2345 Number *values_dofs,
2346 FEEvaluationData<dim, Number, true> &fe_eval,
2347 const bool sum_into_values)
2348 {
2349 const auto &shape_info = fe_eval.get_shape_info();
2350
2351 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
2352 return integrate_tensor_none(n_components,
2353 integration_flag,
2354 values_dofs,
2355 fe_eval,
2356 sum_into_values);
2357 else
2358 return integrate_tensor<fe_degree, n_q_points_1d>(n_components,
2359 integration_flag,
2360 values_dofs,
2361 fe_eval,
2362 sum_into_values);
2363 }
2364 };
2365
2366
2367
2368 template <int dim, typename Number>
2369 struct FEFaceEvaluationImplCollectFromFaceSelector
2370 {
2371 template <int fe_degree>
2372 static bool
2373 run(const unsigned int n_components,
2374 const EvaluationFlags::EvaluationFlags integration_flag,
2375 Number *values_dofs,
2376 FEEvaluationData<dim, Number, true> &fe_eval,
2377 const bool sum_into_values)
2378 {
2379 const auto &shape_info = fe_eval.get_shape_info();
2380 const auto &shape_data = shape_info.data.front();
2381
2382 const unsigned int dofs_per_comp_face =
2383 fe_degree > -1 ?
2384 Utilities::pow(fe_degree + 1, dim - 1) :
2385 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2386
2387 Number *temp = fe_eval.get_scratch_data().begin();
2388 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2389
2390 bool use_vectorization = true;
2391
2392 if (fe_eval.get_dof_access_index() ==
2393 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
2394 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2395 use_vectorization =
2396 fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
2397 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2398 fe_eval.get_cell_ids().end(),
2399 [&](const auto &v) {
2400 return v == fe_eval.get_cell_ids()[0] ||
2401 v == numbers::invalid_unsigned_int;
2402 });
2403
2404 FEFaceEvaluationImplIntegrateSelector<dim, Number>::
2405 template collect_from_face<fe_degree>(n_components,
2406 integration_flag,
2407 values_dofs,
2408 fe_eval,
2409 use_vectorization,
2410 temp,
2411 scratch_data,
2412 sum_into_values);
2413
2414 return false;
2415 }
2416 };
2417
2418
2419
2420 template <int dim, typename Number>
2421 struct FEFaceEvaluationImplIntegrateInFaceSelector
2422 {
2423 template <int fe_degree, int n_q_points_1d>
2424 static bool
2425 run(const unsigned int n_components,
2426 const EvaluationFlags::EvaluationFlags integration_flag,
2427
2428 FEEvaluationData<dim, Number, true> &fe_eval)
2429 {
2430 const auto &shape_info = fe_eval.get_shape_info();
2431 const auto &shape_data = shape_info.data.front();
2432
2433 const unsigned int dofs_per_comp_face =
2434 fe_degree > -1 ?
2435 Utilities::pow(fe_degree + 1, dim - 1) :
2436 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2437
2438 Number *temp = fe_eval.get_scratch_data().begin();
2439 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2440
2441 FEFaceEvaluationImplIntegrateSelector<dim, Number>::
2442 template integrate_in_face<fe_degree, n_q_points_1d>(
2443 n_components, integration_flag, fe_eval, temp, scratch_data);
2444
2445 return false;
2446 }
2447 };
2448
2449
2450
2451 template <int n_face_orientations,
2452 typename Processor,
2453 typename EvaluationData,
2454 const bool check_face_orientations = false>
2455 void
2456 fe_face_evaluation_process_and_io(
2457 Processor &proc,
2458 const unsigned int n_components,
2459 const EvaluationFlags::EvaluationFlags evaluation_flag,
2460 typename Processor::Number2_ *global_vector_ptr,
2461 const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
2462 const EvaluationData &fe_eval,
2463 typename Processor::VectorizedArrayType_ *temp1)
2464 {
2465 constexpr int dim = Processor::dim_;
2466 constexpr int fe_degree = Processor::fe_degree_;
2467 using VectorizedArrayType = typename Processor::VectorizedArrayType_;
2468 constexpr int n_lanes = VectorizedArrayType::size();
2469
2470 using Number = typename Processor::Number_;
2471 using Number2_ = typename Processor::Number2_;
2472
2473 const auto &shape_data = fe_eval.get_shape_info().data.front();
2474 constexpr bool integrate = Processor::do_integrate;
2475 const unsigned int face_no = fe_eval.get_face_no();
2476 const auto &dof_info = fe_eval.get_dof_info();
2477 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
2478 const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
2479 fe_eval.get_dof_access_index();
2480 AssertIndexRange(cell,
2481 dof_info.index_storage_variants[dof_access_index].size());
2482 constexpr unsigned int dofs_per_face =
2483 Utilities::pow(fe_degree + 1, dim - 1);
2484 const unsigned int subface_index = fe_eval.get_subface_index();
2485
2486 const unsigned int n_filled_lanes =
2487 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2488
2489 bool all_faces_are_same = n_filled_lanes == n_lanes;
2490 if (n_face_orientations == n_lanes)
2491 for (unsigned int v = 1; v < n_lanes; ++v)
2492 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
2493 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
2494 {
2495 all_faces_are_same = false;
2496 break;
2497 }
2498
2499 // check for re-orientation ...
2500 std::array<const unsigned int *, n_face_orientations> orientation = {};
2501
2502 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
2503 fe_eval.is_interior_face() == 0)
2504 for (unsigned int v = 0; v < n_lanes; ++v)
2505 {
2506 // the loop breaks once an invalid_unsigned_int is hit for
2507 // all cases except the exterior faces in the ECL loop (where
2508 // some faces might be at the boundaries but others not)
2509 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2510 continue;
2511
2512 if (shape_data.nodal_at_cell_boundaries &&
2513 fe_eval.get_face_orientation(v) != 0)
2514 {
2515 // ... and in case we detect a re-orientation, go to the other
2516 // version of this function that actually allows for this
2517 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2518 check_face_orientations == false)
2519 {
2520 fe_face_evaluation_process_and_io<n_face_orientations,
2521 Processor,
2522 EvaluationData,
2523 true>(proc,
2524 n_components,
2525 evaluation_flag,
2526 global_vector_ptr,
2527 sm_ptr,
2528 fe_eval,
2529 temp1);
2530 return;
2531 }
2532 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2533 fe_eval.get_face_orientation(v), 0);
2534 }
2535 }
2536 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
2537 {
2538 // go to the other version of this function
2539 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2540 check_face_orientations == false)
2541 {
2542 fe_face_evaluation_process_and_io<n_face_orientations,
2543 Processor,
2544 EvaluationData,
2545 true>(proc,
2546 n_components,
2547 evaluation_flag,
2548 global_vector_ptr,
2549 sm_ptr,
2550 fe_eval,
2551 temp1);
2552 return;
2553 }
2554 for (unsigned int v = 0; v < n_face_orientations; ++v)
2555 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2556 fe_eval.get_face_orientation(), 0);
2557 }
2558
2559 // we know that the gradient weights for the Hermite case on the
2560 // right (side==1) are the negative from the value at the left
2561 // (side==0), so we only read out one of them.
2562 VectorizedArrayType grad_weight =
2563 shape_data
2564 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
2565 (1 + face_no % 2))];
2566
2567 // face_to_cell_index_hermite
2568 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2569 {};
2570 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2571 {
2572 if (n_face_orientations == 1)
2573 index_array_hermite[0] =
2574 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
2575 else
2576 {
2577 for (unsigned int v = 0; v < n_lanes; ++v)
2578 {
2579 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2580 continue;
2581
2582 const auto face_no = fe_eval.get_face_no(v);
2583
2584 grad_weight[v] =
2585 shape_data
2586 .shape_data_on_face[0][fe_degree + (integrate ?
2587 (2 - (face_no % 2)) :
2588 (1 + (face_no % 2)))];
2589
2590 index_array_hermite[v] =
2591 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
2592 0);
2593 }
2594 }
2595 }
2596
2597 // face_to_cell_index_nodal
2598 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
2599 {};
2600 if (shape_data.nodal_at_cell_boundaries == true)
2601 {
2602 if (n_face_orientations == 1)
2603 index_array_nodal[0] =
2604 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
2605 else
2606 {
2607 for (unsigned int v = 0; v < n_lanes; ++v)
2608 {
2609 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2610 continue;
2611
2612 const auto face_no = fe_eval.get_face_no(v);
2613
2614 index_array_nodal[v] =
2615 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
2616 0);
2617 }
2618 }
2619 }
2620
2621
2622 const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2623 return (!check_face_orientations || orientation[v] == nullptr) ?
2624 i :
2625 orientation[v][i];
2626 };
2627
2628 const unsigned int cell_index =
2629 dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
2630 fe_eval.get_cell_ids()[0] :
2631 cell * n_lanes;
2632 const unsigned int *dof_indices =
2633 &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
2634
2635 for (unsigned int comp = 0; comp < n_components; ++comp)
2636 {
2637 const std::size_t index_offset =
2638 dof_info.component_dof_indices_offset
2639 [fe_eval.get_active_fe_index()]
2640 [fe_eval.get_first_selected_component()] +
2641 comp * Utilities::pow(fe_degree + 1, dim);
2642
2643 // case 1: contiguous and interleaved indices
2644 if (n_face_orientations == 1 &&
2645 dof_info.index_storage_variants[dof_access_index][cell] ==
2646 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2647 interleaved_contiguous)
2648 {
2649 AssertDimension(
2650 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2651 n_lanes);
2652 Number2_ *vector_ptr =
2653 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
2654
2655 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2656 {
2657 for (unsigned int i = 0; i < dofs_per_face; ++i)
2658 {
2659 Assert(n_face_orientations == 1, ExcNotImplemented());
2660
2661 const unsigned int ind1 = index_array_hermite[0][2 * i];
2662 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
2663 const unsigned int i_ = reorientate(0, i);
2664 proc.hermite_grad_vectorized(temp1[i_],
2665 temp1[i_ + dofs_per_face],
2666 vector_ptr + ind1 * n_lanes,
2667 vector_ptr + ind2 * n_lanes,
2668 grad_weight);
2669 }
2670 }
2671 else
2672 {
2673 for (unsigned int i = 0; i < dofs_per_face; ++i)
2674 {
2675 Assert(n_face_orientations == 1, ExcNotImplemented());
2676
2677 const unsigned int i_ = reorientate(0, i);
2678 const unsigned int ind = index_array_nodal[0][i];
2679 proc.value_vectorized(temp1[i_],
2680 vector_ptr + ind * n_lanes);
2681 }
2682 }
2683 }
2684
2685 // case 2: contiguous and interleaved indices with fixed stride
2686 else if (n_face_orientations == 1 &&
2687 dof_info.index_storage_variants[dof_access_index][cell] ==
2688 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2689 interleaved_contiguous_strided)
2690 {
2691 AssertDimension(
2692 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2693 n_lanes);
2694 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
2695 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2696 {
2697 for (unsigned int i = 0; i < dofs_per_face; ++i)
2698 {
2699 Assert(n_face_orientations == 1, ExcNotImplemented());
2700
2701 const unsigned int i_ = reorientate(0, i);
2702 const unsigned int ind1 =
2703 index_array_hermite[0][2 * i] * n_lanes;
2704 const unsigned int ind2 =
2705 index_array_hermite[0][2 * i + 1] * n_lanes;
2706 proc.hermite_grad_vectorized_indexed(
2707 temp1[i_],
2708 temp1[i_ + dofs_per_face],
2709 vector_ptr + ind1,
2710 vector_ptr + ind2,
2711 grad_weight,
2712 dof_indices,
2713 dof_indices);
2714 }
2715 }
2716 else
2717 {
2718 for (unsigned int i = 0; i < dofs_per_face; ++i)
2719 {
2720 Assert(n_face_orientations == 1, ExcNotImplemented());
2721
2722 const unsigned int i_ = reorientate(0, i);
2723 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
2724 proc.value_vectorized_indexed(temp1[i_],
2725 vector_ptr + ind,
2726 dof_indices);
2727 }
2728 }
2729 }
2730
2731 // case 3: contiguous and interleaved indices with mixed stride
2732 else if (n_face_orientations == 1 &&
2733 dof_info.index_storage_variants[dof_access_index][cell] ==
2734 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2735 interleaved_contiguous_mixed_strides)
2736 {
2737 const unsigned int *strides =
2738 &dof_info.dof_indices_interleave_strides[dof_access_index]
2739 [cell * n_lanes];
2740 unsigned int indices[n_lanes];
2741 for (unsigned int v = 0; v < n_lanes; ++v)
2742 indices[v] = dof_indices[v] + index_offset * strides[v];
2743 const unsigned int n_filled_lanes =
2744 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2745
2746 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2747 {
2748 if (n_filled_lanes == n_lanes)
2749 for (unsigned int i = 0; i < dofs_per_face; ++i)
2750 {
2751 Assert(n_face_orientations == 1, ExcNotImplemented());
2752
2753 const unsigned int i_ = reorientate(0, i);
2754 unsigned int ind1[n_lanes];
2755 DEAL_II_OPENMP_SIMD_PRAGMA
2756 for (unsigned int v = 0; v < n_lanes; ++v)
2757 ind1[v] = indices[v] +
2758 index_array_hermite[0][2 * i] * strides[v];
2759 unsigned int ind2[n_lanes];
2760 DEAL_II_OPENMP_SIMD_PRAGMA
2761 for (unsigned int v = 0; v < n_lanes; ++v)
2762 ind2[v] =
2763 indices[v] +
2764 // TODO
2765 index_array_hermite[0][2 * i + 1] * strides[v];
2766 proc.hermite_grad_vectorized_indexed(
2767 temp1[i_],
2768 temp1[i_ + dofs_per_face],
2769 global_vector_ptr,
2770 global_vector_ptr,
2771 grad_weight,
2772 ind1,
2773 ind2);
2774 }
2775 else
2776 {
2777 if (integrate == false)
2778 for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
2779 temp1[i] = VectorizedArrayType();
2780
2781 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2782 for (unsigned int i = 0; i < dofs_per_face; ++i)
2783 {
2784 const unsigned int i_ =
2785 reorientate(n_face_orientations == 1 ? 0 : v, i);
2786 proc.hermite_grad(
2787 temp1[i_][v],
2788 temp1[i_ + dofs_per_face][v],
2789 global_vector_ptr
2790 [indices[v] +
2791 index_array_hermite
2792 [n_face_orientations == 1 ? 0 : v][2 * i] *
2793 strides[v]],
2794 global_vector_ptr
2795 [indices[v] +
2796 index_array_hermite[n_face_orientations == 1 ?
2797 0 :
2798 v][2 * i + 1] *
2799 strides[v]],
2800 grad_weight[n_face_orientations == 1 ? 0 : v]);
2801 }
2802 }
2803 }
2804 else
2805 {
2806 if (n_filled_lanes == n_lanes)
2807 for (unsigned int i = 0; i < dofs_per_face; ++i)
2808 {
2809 Assert(n_face_orientations == 1, ExcInternalError());
2810 unsigned int ind[n_lanes];
2811 DEAL_II_OPENMP_SIMD_PRAGMA
2812 for (unsigned int v = 0; v < n_lanes; ++v)
2813 ind[v] =
2814 indices[v] + index_array_nodal[0][i] * strides[v];
2815 const unsigned int i_ = reorientate(0, i);
2816 proc.value_vectorized_indexed(temp1[i_],
2817 global_vector_ptr,
2818 ind);
2819 }
2820 else
2821 {
2822 if (integrate == false)
2823 for (unsigned int i = 0; i < dofs_per_face; ++i)
2824 temp1[i] = VectorizedArrayType();
2825
2826 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2827 for (unsigned int i = 0; i < dofs_per_face; ++i)
2828 proc.value(
2829 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
2830 i)][v],
2831 global_vector_ptr
2832 [indices[v] +
2833 index_array_nodal[n_face_orientations == 1 ? 0 : v]
2834 [i] *
2835 strides[v]]);
2836 }
2837 }
2838 }
2839
2840 // case 4: contiguous indices without interleaving
2841 else if (n_face_orientations > 1 ||
2842 dof_info.index_storage_variants[dof_access_index][cell] ==
2843 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2844 contiguous)
2845 {
2846 Number2_ *vector_ptr = global_vector_ptr + index_offset;
2847
2848 const bool vectorization_possible =
2849 all_faces_are_same && (sm_ptr == nullptr);
2850
2851 std::array<Number2_ *, n_lanes> vector_ptrs{{nullptr}};
2852 std::array<unsigned int, n_lanes> reordered_indices{
2853 {numbers::invalid_unsigned_int}};
2854
2855 if (vectorization_possible == false)
2856 {
2857 if (n_face_orientations == 1)
2858 {
2859 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2860 if (sm_ptr == nullptr)
2861 {
2862 vector_ptrs[v] = vector_ptr + dof_indices[v];
2863 }
2864 else
2865 {
2866 const auto &temp =
2867 dof_info
2868 .dof_indices_contiguous_sm[dof_access_index]
2869 [cell * n_lanes + v];
2870 vector_ptrs[v] = const_cast<Number2_ *>(
2871 sm_ptr->operator[](temp.first).data() +
2872 temp.second + index_offset);
2873 }
2874 }
2875 else if (n_face_orientations == n_lanes)
2876 {
2877 const auto &cells = fe_eval.get_cell_ids();
2878 for (unsigned int v = 0; v < n_lanes; ++v)
2879 if (cells[v] != numbers::invalid_unsigned_int)
2880 {
2881 if (sm_ptr == nullptr)
2882 {
2883 vector_ptrs[v] =
2884 vector_ptr +
2885 dof_info
2886 .dof_indices_contiguous[dof_access_index]
2887 [cells[v]];
2888 }
2889 else
2890 {
2891 const auto &temp =
2892 dof_info
2893 .dof_indices_contiguous_sm[dof_access_index]
2894 [cells[v]];
2895 vector_ptrs[v] = const_cast<Number2_ *>(
2896 sm_ptr->operator[](temp.first).data() +
2897 temp.second + index_offset);
2898 }
2899 }
2900 }
2901 else
2902 {
2903 DEAL_II_NOT_IMPLEMENTED();
2904 }
2905 }
2906 else if (n_face_orientations == n_lanes)
2907 {
2908 for (unsigned int v = 0; v < n_lanes; ++v)
2909 reordered_indices[v] =
2910 dof_info.dof_indices_contiguous[dof_access_index]
2911 [fe_eval.get_cell_ids()[v]];
2912 dof_indices = reordered_indices.data();
2913 }
2914
2915 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2916 {
2917 if (vectorization_possible)
2918 for (unsigned int i = 0; i < dofs_per_face; ++i)
2919 {
2920 const unsigned int ind1 = index_array_hermite[0][2 * i];
2921 const unsigned int ind2 =
2922 index_array_hermite[0][2 * i + 1];
2923 const unsigned int i_ = reorientate(0, i);
2924
2925 proc.hermite_grad_vectorized_indexed(
2926 temp1[i_],
2927 temp1[i_ + dofs_per_face],
2928 vector_ptr + ind1,
2929 vector_ptr + ind2,
2930 grad_weight,
2931 dof_indices,
2932 dof_indices);
2933 }
2934 else if (n_face_orientations == 1)
2935 for (unsigned int i = 0; i < dofs_per_face; ++i)
2936 {
2937 const unsigned int ind1 = index_array_hermite[0][2 * i];
2938 const unsigned int ind2 =
2939 index_array_hermite[0][2 * i + 1];
2940 const unsigned int i_ = reorientate(0, i);
2941
2942 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2943 proc.hermite_grad(temp1[i_][v],
2944 temp1[i_ + dofs_per_face][v],
2945 vector_ptrs[v][ind1],
2946 vector_ptrs[v][ind2],
2947 grad_weight[v]);
2948
2949 if (integrate == false)
2950 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2951 {
2952 temp1[i][v] = 0.0;
2953 temp1[i + dofs_per_face][v] = 0.0;
2954 }
2955 }
2956 else
2957 {
2958 if (integrate == false && n_filled_lanes < n_lanes)
2959 for (unsigned int i = 0; i < dofs_per_face; ++i)
2960 temp1[i] = temp1[i + dofs_per_face] = Number();
2961
2962 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2963 for (unsigned int i = 0; i < dofs_per_face; ++i)
2964 proc.hermite_grad(
2965 temp1[reorientate(v, i)][v],
2966 temp1[reorientate(v, i) + dofs_per_face][v],
2967 vector_ptrs[v][index_array_hermite[v][2 * i]],
2968 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
2969 grad_weight[v]);
2970 }
2971 }
2972 else
2973 {
2974 if (vectorization_possible)
2975 for (unsigned int i = 0; i < dofs_per_face; ++i)
2976 {
2977 const unsigned int ind = index_array_nodal[0][i];
2978 const unsigned int i_ = reorientate(0, i);
2979
2980 proc.value_vectorized_indexed(temp1[i_],
2981 vector_ptr + ind,
2982 dof_indices);
2983 }
2984 else
2985 {
2986 if constexpr (n_face_orientations == 1)
2987 for (unsigned int i = 0; i < dofs_per_face; ++i)
2988 {
2989 const unsigned int ind = index_array_nodal[0][i];
2990 const unsigned int i_ = reorientate(0, i);
2991
2992 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2993 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
2994
2995 if constexpr (integrate == false)
2996 for (unsigned int v = n_filled_lanes; v < n_lanes;
2997 ++v)
2998 temp1[i_][v] = 0.0;
2999 }
3000 else
3001 {
3002 if (integrate == false && n_filled_lanes < n_lanes)
3003 for (unsigned int i = 0; i < dofs_per_face; ++i)
3004 temp1[i] = Number();
3005
3006 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3007 for (unsigned int i = 0; i < dofs_per_face; ++i)
3008 proc.value(temp1[reorientate(v, i)][v],
3009 vector_ptrs[v][index_array_nodal[v][i]]);
3010 }
3011 }
3012 }
3013 }
3014 else
3015 {
3016 // We should not end up here, this should be caught by
3017 // FEFaceEvaluationImplGatherEvaluateSelector::supports()
3018 DEAL_II_ASSERT_UNREACHABLE();
3019 }
3020 temp1 += 3 * dofs_per_face;
3021 }
3022 }
3023
3024
3025
3026 template <int dim, typename Number2, typename VectorizedArrayType>
3027 struct FEFaceEvaluationImplGatherEvaluateSelector
3028 {
3029 using Number = typename VectorizedArrayType::value_type;
3030
3031 template <int fe_degree, int n_q_points_1d>
3032 static bool
3033 run(const unsigned int n_components,
3034 const EvaluationFlags::EvaluationFlags evaluation_flag,
3035 const Number2 *src_ptr,
3036 const std::vector<ArrayView<const Number2>> *sm_ptr,
3037 FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
3038 {
3039 Assert(fe_degree > -1, ExcInternalError());
3040 Assert(fe_eval.get_shape_info().element_type <=
3041 MatrixFreeFunctions::tensor_symmetric_no_collocation,
3042 ExcInternalError());
3043
3044 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3045
3046 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3047 VectorizedArrayType *scratch_data =
3048 temp + 3 * n_components * dofs_per_face;
3049
3050 Processor<fe_degree> p;
3051
3052 if (fe_eval.get_dof_access_index() ==
3053 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3054 fe_eval.is_interior_face() == false)
3055 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3056 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
3057 else
3058 fe_face_evaluation_process_and_io<1>(
3059 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
3060
3061 const unsigned int subface_index = fe_eval.get_subface_index();
3062
3063 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3064 FEFaceEvaluationImpl<true,
3065 dim,
3066 fe_degree,
3067 n_q_points_1d,
3068 VectorizedArrayType>::
3069 evaluate_in_face(n_components,
3070 evaluation_flag,
3071 fe_eval.get_shape_info().data.front(),
3072 temp,
3073 fe_eval.begin_values(),
3074 fe_eval.begin_gradients(),
3075 fe_eval.begin_hessians(),
3076 scratch_data,
3077 subface_index);
3078 else
3079 FEFaceEvaluationImpl<false,
3080 dim,
3081 fe_degree,
3082 n_q_points_1d,
3083 VectorizedArrayType>::
3084 evaluate_in_face(n_components,
3085 evaluation_flag,
3086 fe_eval.get_shape_info().data.front(),
3087 temp,
3088 fe_eval.begin_values(),
3089 fe_eval.begin_gradients(),
3090 fe_eval.begin_hessians(),
3091 scratch_data,
3092 subface_index);
3093
3094 // re-orientation for cases not possible with above algorithm
3095 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3096 {
3097 if (fe_eval.get_dof_access_index() ==
3098 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3099 fe_eval.is_interior_face() == false)
3100 {
3101 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3102 {
3103 // the loop breaks once an invalid_unsigned_int is hit for
3104 // all cases except the exterior faces in the ECL loop (where
3105 // some faces might be at the boundaries but others not)
3106 if (fe_eval.get_cell_ids()[v] ==
3107 numbers::invalid_unsigned_int)
3108 continue;
3109
3110 if (fe_eval.get_face_orientation(v) != 0)
3111 adjust_for_face_orientation_per_lane(
3112 dim,
3113 n_components,
3114 v,
3115 evaluation_flag,
3116 &fe_eval.get_shape_info().face_orientations_quad(
3117 fe_eval.get_face_orientation(v), 0),
3118 false,
3119 Utilities::pow(n_q_points_1d, dim - 1),
3120 &temp[0][0],
3121 fe_eval.begin_values(),
3122 fe_eval.begin_gradients(),
3123 fe_eval.begin_hessians());
3124 }
3125 }
3126 else if (fe_eval.get_face_orientation() != 0)
3127 adjust_for_face_orientation(
3128 dim,
3129 n_components,
3130 evaluation_flag,
3131 &fe_eval.get_shape_info().face_orientations_quad(
3132 fe_eval.get_face_orientation(), 0),
3133 false,
3134 Utilities::pow(n_q_points_1d, dim - 1),
3135 temp,
3136 fe_eval.begin_values(),
3137 fe_eval.begin_gradients(),
3138 fe_eval.begin_hessians());
3139 }
3140
3141 return false;
3142 }
3143
3144 template <typename Number3>
3145 static bool
3146 supports(const EvaluationFlags::EvaluationFlags evaluation_flag,
3147 const MatrixFreeFunctions::ShapeInfo<Number3> &shape_info,
3148 const Number2 *vector_ptr,
3149 MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
3150 {
3151 const unsigned int fe_degree = shape_info.data.front().fe_degree;
3152 if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
3153 (evaluation_flag & EvaluationFlags::gradients &&
3154 (fe_degree < 2 ||
3155 shape_info.data.front().element_type !=
3156 MatrixFreeFunctions::tensor_symmetric_hermite)) ||
3157 (evaluation_flag & EvaluationFlags::hessians) ||
3158 vector_ptr == nullptr ||
3159 shape_info.data.front().element_type >
3160 MatrixFreeFunctions::tensor_symmetric_no_collocation ||
3161 storage <
3162 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
3163 return false;
3164 else
3165 return true;
3166 }
3167
3168 private:
3169 template <int fe_degree>
3170 struct Processor
3171 {
3172 static const bool do_integrate = false;
3173 static const int dim_ = dim;
3174 static const int fe_degree_ = fe_degree;
3175 using VectorizedArrayType_ = VectorizedArrayType;
3176 using Number_ = Number;
3177 using Number2_ = const Number2;
3178
3179 template <typename T0, typename T1, typename T2>
3180 void
3181 hermite_grad_vectorized(T0 &temp_1,
3182 T0 &temp_2,
3183 const T1 src_ptr_1,
3184 const T1 src_ptr_2,
3185 const T2 &grad_weight)
3186 {
3187 do_vectorized_read(src_ptr_1, temp_1);
3188 do_vectorized_read(src_ptr_2, temp_2);
3189 temp_2 = grad_weight * (temp_1 - temp_2);
3190 }
3191
3192 template <typename T1, typename T2>
3193 void
3194 value_vectorized(T1 &temp, const T2 src_ptr)
3195 {
3196 do_vectorized_read(src_ptr, temp);
3197 }
3198
3199 template <typename T0, typename T1, typename T2, typename T3>
3200 void
3201 hermite_grad_vectorized_indexed(T0 &temp_1,
3202 T0 &temp_2,
3203 const T1 src_ptr_1,
3204 const T1 src_ptr_2,
3205 const T2 &grad_weight,
3206 const T3 &indices_1,
3207 const T3 &indices_2)
3208 {
3209 do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3210 do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3211 temp_2 = grad_weight * (temp_1 - temp_2);
3212 }
3213
3214 template <typename T0, typename T1, typename T2>
3215 void
3216 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
3217 {
3218 do_vectorized_gather(src_ptr, indices, temp);
3219 }
3220
3221 template <typename T0, typename T1, typename T2>
3222 void
3223 hermite_grad(T0 &temp_1,
3224 T0 &temp_2,
3225 const T1 &src_ptr_1,
3226 const T1 &src_ptr_2,
3227 const T2 &grad_weight)
3228 {
3229 // case 3a)
3230 temp_1 = src_ptr_1;
3231 temp_2 = grad_weight * (temp_1 - src_ptr_2);
3232 }
3233
3234 template <typename T1, typename T2>
3235 void
3236 value(T1 &temp, const T2 &src_ptr)
3237 {
3238 // case 3b)
3239 temp = src_ptr;
3240 }
3241 };
3242 };
3243
3244
3245
3246 template <int dim, typename Number2, typename VectorizedArrayType>
3247 struct FEFaceEvaluationImplIntegrateScatterSelector
3248 {
3249 using Number = typename VectorizedArrayType::value_type;
3250
3251 template <int fe_degree, int n_q_points_1d>
3252 static bool
3253 run(const unsigned int n_components,
3254 const EvaluationFlags::EvaluationFlags integration_flag,
3255 Number2 *dst_ptr,
3256 const std::vector<ArrayView<const Number2>> *sm_ptr,
3257 FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
3258 {
3259 Assert(fe_degree > -1, ExcInternalError());
3260 Assert(fe_eval.get_shape_info().element_type <=
3261 MatrixFreeFunctions::tensor_symmetric_no_collocation,
3262 ExcInternalError());
3263
3264 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3265
3266 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3267 VectorizedArrayType *scratch_data =
3268 temp + 3 * n_components * dofs_per_face;
3269
3270 const unsigned int subface_index = fe_eval.get_subface_index();
3271
3272 // re-orientation for cases not possible with the io function below
3273 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3274 {
3275 if (fe_eval.get_dof_access_index() ==
3276 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3277 fe_eval.is_interior_face() == false)
3278 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3279 {
3280 // the loop breaks once an invalid_unsigned_int is hit for
3281 // all cases except the exterior faces in the ECL loop (where
3282 // some faces might be at the boundaries but others not)
3283 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
3284 continue;
3285
3286 if (fe_eval.get_face_orientation(v) != 0)
3287 adjust_for_face_orientation_per_lane(
3288 dim,
3289 n_components,
3290 v,
3291 integration_flag,
3292 &fe_eval.get_shape_info().face_orientations_quad(
3293 fe_eval.get_face_orientation(v), 0),
3294 true,
3295 Utilities::pow(n_q_points_1d, dim - 1),
3296 &temp[0][0],
3297 fe_eval.begin_values(),
3298 fe_eval.begin_gradients(),
3299 fe_eval.begin_hessians());
3300 }
3301 else if (fe_eval.get_face_orientation() != 0)
3302 adjust_for_face_orientation(
3303 dim,
3304 n_components,
3305 integration_flag,
3306 &fe_eval.get_shape_info().face_orientations_quad(
3307 fe_eval.get_face_orientation(), 0),
3308 true,
3309 Utilities::pow(n_q_points_1d, dim - 1),
3310 temp,
3311 fe_eval.begin_values(),
3312 fe_eval.begin_gradients(),
3313 fe_eval.begin_hessians());
3314 }
3315
3316 if (fe_degree > -1 && fe_eval.get_subface_index() >=
3317 GeometryInfo<dim - 1>::max_children_per_cell)
3318 FEFaceEvaluationImpl<true,
3319 dim,
3320 fe_degree,
3321 n_q_points_1d,
3322 VectorizedArrayType>::
3323 integrate_in_face(n_components,
3324 integration_flag,
3325 fe_eval.get_shape_info().data.front(),
3326 temp,
3327 fe_eval.begin_values(),
3328 fe_eval.begin_gradients(),
3329 fe_eval.begin_hessians(),
3330 scratch_data,
3331 subface_index);
3332 else
3333 FEFaceEvaluationImpl<false,
3334 dim,
3335 fe_degree,
3336 n_q_points_1d,
3337 VectorizedArrayType>::
3338 integrate_in_face(n_components,
3339 integration_flag,
3340 fe_eval.get_shape_info().data.front(),
3341 temp,
3342 fe_eval.begin_values(),
3343 fe_eval.begin_gradients(),
3344 fe_eval.begin_hessians(),
3345 scratch_data,
3346 subface_index);
3347
3348 Processor<fe_degree> p;
3349
3350 if (fe_eval.get_dof_access_index() ==
3351 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3352 fe_eval.is_interior_face() == false)
3353 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3354 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3355 else
3356 fe_face_evaluation_process_and_io<1>(
3357 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3358
3359 return false;
3360 }
3361
3362 private:
3363 template <int fe_degree>
3364 struct Processor
3365 {
3366 static const bool do_integrate = true;
3367 static const int dim_ = dim;
3368 static const int fe_degree_ = fe_degree;
3369 using VectorizedArrayType_ = VectorizedArrayType;
3370 using Number_ = Number;
3371 using Number2_ = Number2;
3372
3373 template <typename T0, typename T1, typename T2, typename T3, typename T4>
3374 void
3375 hermite_grad_vectorized(const T0 &temp_1,
3376 const T1 &temp_2,
3377 T2 dst_ptr_1,
3378 T3 dst_ptr_2,
3379 const T4 &grad_weight)
3380 {
3381 // case 1a)
3382 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3383 const VectorizedArrayType grad = grad_weight * temp_2;
3384 do_vectorized_add(val, dst_ptr_1);
3385 do_vectorized_add(grad, dst_ptr_2);
3386 }
3387
3388 template <typename T0, typename T1>
3389 void
3390 value_vectorized(const T0 &temp, T1 dst_ptr)
3391 {
3392 // case 1b)
3393 do_vectorized_add(temp, dst_ptr);
3394 }
3395
3396 template <typename T0, typename T1, typename T2, typename T3>
3397 void
3398 hermite_grad_vectorized_indexed(const T0 &temp_1,
3399 const T0 &temp_2,
3400 T1 dst_ptr_1,
3401 T1 dst_ptr_2,
3402 const T2 &grad_weight,
3403 const T3 &indices_1,
3404 const T3 &indices_2)
3405 {
3406 // case 2a)
3407 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3408 const VectorizedArrayType grad = grad_weight * temp_2;
3409 do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
3410 do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
3411 }
3412
3413 template <typename T0, typename T1, typename T2>
3414 void
3415 value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
3416 {
3417 // case 2b)
3418 do_vectorized_scatter_add(temp, indices, dst_ptr);
3419 }
3420
3421 template <typename T0, typename T1, typename T2>
3422 void
3423 hermite_grad(const T0 &temp_1,
3424 const T0 &temp_2,
3425 T1 &dst_ptr_1,
3426 T1 &dst_ptr_2,
3427 const T2 &grad_weight)
3428 {
3429 // case 3a)
3430 const Number val = temp_1 - grad_weight * temp_2;
3431 const Number grad = grad_weight * temp_2;
3432 dst_ptr_1 += val;
3433 dst_ptr_2 += grad;
3434 }
3435
3436 template <typename T0, typename T1>
3437 void
3438 value(const T0 &temp, T1 &dst_ptr)
3439 {
3440 // case 3b)
3441 dst_ptr += temp;
3442 }
3443 };
3444 };
3445} // end of namespace internal
3446
3447
3448DEAL_II_NAMESPACE_CLOSE
3449
3450#endif
ArrayView::begin
iterator begin() const
Definition array_view.h:702
FEEvaluationData::get_face_no
std::uint8_t get_face_no(const unsigned int v=0) const
FEEvaluationData::get_dof_access_index
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
FEEvaluationData::shape_info_number_type
ScalarNumber shape_info_number_type
Definition fe_evaluation_data.h:127
FEEvaluationData::get_shape_info
const ShapeInfoType & get_shape_info() const
FEEvaluationData::get_cell_ids
const std::array< unsigned int, n_lanes > & get_cell_ids() const
Definition fe_evaluation_data.h:496
FEEvaluationData::begin_gradients
const Number * begin_gradients() const
FEEvaluationData::get_subface_index
unsigned int get_subface_index() const
FEEvaluationData::is_interior_face
bool is_interior_face() const
FEEvaluationData::get_scratch_data
ArrayView< Number > get_scratch_data() const
FEEvaluationData::begin_values
const Number * begin_values() const
FEEvaluationData::get_face_orientation
std::uint8_t get_face_orientation(const unsigned int v=0) const
FEEvaluationData::begin_hessians
const Number * begin_hessians() const
VectorizedArray::gather
void gather(const Number *base_ptr, const unsigned int *offsets)
Definition vectorization.h:672
VectorizedArray::load
void load(const OtherNumber *ptr)
Definition vectorization.h:587
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:141
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
evaluation_flags.h
DoFCouplingType::contiguous
@ contiguous
fe_evaluation_data.h
geometry_info.h
cell_index
unsigned int cell_index
Definition grid_tools_topology.cc:783
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1780
internal::MatrixFreeFunctions::tensor_symmetric_no_collocation
@ tensor_symmetric_no_collocation
Definition shape_info.h:87
internal::MatrixFreeFunctions::tensor_raviart_thomas
@ tensor_raviart_thomas
Definition shape_info.h:112
internal::MatrixFreeFunctions::tensor_none
@ tensor_none
Definition shape_info.h:117
internal::MatrixFreeFunctions::tensor_symmetric
@ tensor_symmetric
Definition shape_info.h:81
internal::MatrixFreeFunctions::tensor_symmetric_hermite
@ tensor_symmetric_hermite
Definition shape_info.h:74
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exceptions.h:1938
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
dof_info.h
data
std::vector< index_type > data
Definition mpi.cc:735
EvaluationFlags::EvaluationFlags
EvaluationFlags
The EvaluationFlags enum.
Definition evaluation_flags.h:42
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::do_vectorized_add
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
Definition evaluation_kernels_face.h:1260
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition tensor_product_kernels.h:2433
internal::adjust_for_face_orientation_per_lane
void adjust_for_face_orientation_per_lane(const unsigned int dim, const unsigned int n_components, const unsigned int v, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad=nullptr, VectorizedArrayType *hessians_quad=nullptr)
Definition evaluation_kernels_face.h:1379
internal::fe_face_evaluation_process_and_io
void fe_face_evaluation_process_and_io(Processor &proc, const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, typename Processor::Number2_ *global_vector_ptr, const std::vector< ArrayView< const typename Processor::Number2_ > > *sm_ptr, const EvaluationData &fe_eval, typename Processor::VectorizedArrayType_ *temp1)
Definition evaluation_kernels_face.h:2456
internal::apply_matrix_vector_product
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
Definition tensor_product_kernels.h:118
internal::do_vectorized_scatter_add
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
Definition evaluation_kernels_face.h:1284
internal::do_vectorized_gather
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
Definition evaluation_kernels_face.h:1234
internal::do_vectorized_read
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
Definition evaluation_kernels_face.h:1212
internal::interpolate_to_face
std::enable_if_t< contract_onto_face, void > interpolate_to_face(const Number2 *shape_values, const std::array< int, 2 > &n_blocks, const std::array< int, 2 > &steps, const Number *input, Number *DEAL_II_RESTRICT output, const int n_rows_runtime=0, const int stride_runtime=1)
Definition tensor_product_kernels.h:2365
internal::adjust_for_face_orientation
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad, Number *hessians_quad)
Definition evaluation_kernels_face.h:1316
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
shape_info.h
internal::FEFaceEvaluationImplCollectFromFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
Definition evaluation_kernels_face.h:2373
internal::FEFaceEvaluationImplEvaluateInFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1873
internal::FEFaceEvaluationImplEvaluateSelector::evaluate_in_face
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval, Number *temp, Number *scratch_data)
Definition evaluation_kernels_face.h:1632
internal::FEFaceEvaluationImplEvaluateSelector::project_to_face
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, Number *temp, Number *scratch_data)
Definition evaluation_kernels_face.h:1572
internal::FEFaceEvaluationImplEvaluateSelector::evaluate_tensor_none
static bool evaluate_tensor_none(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1450
internal::FEFaceEvaluationImplEvaluateSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1799
internal::FEFaceEvaluationImplEvaluateSelector::adjust_quadrature_for_face_orientation
static void adjust_quadrature_for_face_orientation(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, Number *temp)
Definition evaluation_kernels_face.h:1699
internal::FEFaceEvaluationImplEvaluateSelector::evaluate_tensor
static bool evaluate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1752
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::hermite_grad
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T1 &src_ptr_2, const T2 &grad_weight)
Definition evaluation_kernels_face.h:3223
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::value_vectorized_indexed
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
Definition evaluation_kernels_face.h:3216
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::hermite_grad_vectorized
void hermite_grad_vectorized(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight)
Definition evaluation_kernels_face.h:3181
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::hermite_grad_vectorized_indexed
void hermite_grad_vectorized_indexed(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
Definition evaluation_kernels_face.h:3201
internal::FEFaceEvaluationImplGatherEvaluateSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number2 *src_ptr, const std::vector< ArrayView< const Number2 > > *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
Definition evaluation_kernels_face.h:3033
internal::FEFaceEvaluationImplGatherEvaluateSelector::supports
static bool supports(const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number3 > &shape_info, const Number2 *vector_ptr, MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
Definition evaluation_kernels_face.h:3146
internal::FEFaceEvaluationImplIntegrateInFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:2425
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::hermite_grad_vectorized
void hermite_grad_vectorized(const T0 &temp_1, const T1 &temp_2, T2 dst_ptr_1, T3 dst_ptr_2, const T4 &grad_weight)
Definition evaluation_kernels_face.h:3375
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::hermite_grad_vectorized_indexed
void hermite_grad_vectorized_indexed(const T0 &temp_1, const T0 &temp_2, T1 dst_ptr_1, T1 dst_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
Definition evaluation_kernels_face.h:3398
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::hermite_grad
void hermite_grad(const T0 &temp_1, const T0 &temp_2, T1 &dst_ptr_1, T1 &dst_ptr_2, const T2 &grad_weight)
Definition evaluation_kernels_face.h:3423
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::value_vectorized_indexed
void value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
Definition evaluation_kernels_face.h:3415
internal::FEFaceEvaluationImplIntegrateScatterSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number2 *dst_ptr, const std::vector< ArrayView< const Number2 > > *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
Definition evaluation_kernels_face.h:3253
internal::FEFaceEvaluationImplIntegrateSelector::integrate_tensor_none
static bool integrate_tensor_none(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
Definition evaluation_kernels_face.h:1904
internal::FEFaceEvaluationImplIntegrateSelector::integrate_tensor
static bool integrate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
Definition evaluation_kernels_face.h:2291
internal::FEFaceEvaluationImplIntegrateSelector::collect_from_face
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, const Number *temp, Number *scratch_data, const bool sum_into_values)
Definition evaluation_kernels_face.h:2224
internal::FEFaceEvaluationImplIntegrateSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
Definition evaluation_kernels_face.h:2343
internal::FEFaceEvaluationImplIntegrateSelector::integrate_in_face
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval, Number *temp, Number *scratch_data)
Definition evaluation_kernels_face.h:2155
internal::FEFaceEvaluationImplIntegrateSelector::adjust_quadrature_for_face_orientation
static void adjust_quadrature_for_face_orientation(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, Number *temp)
Definition evaluation_kernels_face.h:2100
internal::FEFaceEvaluationImplProjectToFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1826
internal::FEFaceEvaluationImplRaviartThomas::evaluate_or_integrate_in_face
static void evaluate_or_integrate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag, const std::vector< MatrixFreeFunctions::UnivariateShapeData< Number2 > > &shape_data, Number *values_dofs_in, Number *values, Number *gradients, Number *scratch_data, const unsigned int subface_index, const unsigned int face_direction)
Definition evaluation_kernels_face.h:534
internal::FEFaceEvaluationImplRaviartThomas::Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
Definition evaluation_kernels_face.h:526
internal::FEFaceEvaluationImpl::Eval
EvaluatorTensorProduct< symmetric_evaluate ? evaluate_evenodd :evaluate_general, dim - 1, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
Definition evaluation_kernels_face.h:59
internal::FEFaceEvaluationImpl::evaluate_in_face
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
Definition evaluation_kernels_face.h:92
internal::FEFaceEvaluationImpl::create_evaluator_tensor_product
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const unsigned int subface_index, const unsigned int direction)
Definition evaluation_kernels_face.h:62
internal::FEFaceEvaluationImpl::integrate_in_face
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
Definition evaluation_kernels_face.h:290
internal::FEFaceEvaluationImpl::Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
Definition evaluation_kernels_face.h:51
internal::FEFaceNormalEvaluationImpl::Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
Definition evaluation_kernels_face.h:858
internal::FEFaceNormalEvaluationImpl::interpolate_quadrature
static void interpolate_quadrature(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const Number *input, Number *output, const unsigned int face_no)
Definition evaluation_kernels_face.h:894
internal::FEFaceNormalEvaluationImpl::interpolate_generic
static void interpolate_generic(const unsigned int n_components, const Number *input, Number *output, const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const unsigned int n_points_1d, const std::array< AlignedVector< Number2 >, 2 > &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face)
Definition evaluation_kernels_face.h:922
internal::FEFaceNormalEvaluationImpl::interpolate
static void interpolate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const Number *input, Number *output, const unsigned int face_no)
Definition evaluation_kernels_face.h:862
internal::FEFaceNormalEvaluationImpl::interpolate_raviart_thomas
static void interpolate_raviart_thomas(const unsigned int n_components, const Number *input, Number *output, const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info)
Definition evaluation_kernels_face.h:1023
internal::MatrixFreeFunctions::ShapeInfo::data
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:485
internal::MatrixFreeFunctions::ShapeInfo::face_orientations_quad
::Table< 2, unsigned int > face_orientations_quad
Definition shape_info.h:624
tensor_product_kernels.h
Generated by doxygen 1.9.8