Reference documentation for deal.II version 9.4.1
\(\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.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2017 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_evaluation_kernels_h
18#define dealii_matrix_free_evaluation_kernels_h
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/base/geometry_info.h>
23#include <deal.II/base/ndarray.h>
24#include <deal.II/base/utilities.h>
25#include <deal.II/base/vectorization.h>
26
27#include <deal.II/matrix_free/dof_info.h>
28#include <deal.II/matrix_free/evaluation_flags.h>
29#include <deal.II/matrix_free/fe_evaluation_data.h>
30#include <deal.II/matrix_free/shape_info.h>
31#include <deal.II/matrix_free/tensor_product_kernels.h>
32
33
34DEAL_II_NAMESPACE_OPEN
35
36
37namespace internal
38{
39 // Select evaluator type from element shape function type
40 template <MatrixFreeFunctions::ElementType element, bool is_long>
41 struct EvaluatorSelector
42 {};
43
44 template <bool is_long>
45 struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
46 {
47 static const EvaluatorVariant variant = evaluate_general;
48 };
49
50 template <>
51 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
52 {
53 static const EvaluatorVariant variant = evaluate_symmetric;
54 };
55
56 template <>
57 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
58 {
59 static const EvaluatorVariant variant = evaluate_evenodd;
60 };
61
62 template <bool is_long>
63 struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
64 {
65 static const EvaluatorVariant variant = evaluate_general;
66 };
67
68 template <>
69 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
70 false>
71 {
72 static const EvaluatorVariant variant = evaluate_general;
73 };
74
75 template <>
76 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
77 {
78 static const EvaluatorVariant variant = evaluate_evenodd;
79 };
80
81 template <bool is_long>
82 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
83 is_long>
84 {
85 static const EvaluatorVariant variant = evaluate_evenodd;
86 };
87
88 template <bool is_long>
89 struct EvaluatorSelector<MatrixFreeFunctions::tensor_raviart_thomas, is_long>
90 {
91 static const EvaluatorVariant variant = evaluate_raviart_thomas;
92 };
93
94
95
112 template <MatrixFreeFunctions::ElementType type,
113 int dim,
114 int fe_degree,
115 int n_q_points_1d,
116 typename Number>
117 struct FEEvaluationImpl
118 {
119 static const EvaluatorVariant variant =
120 EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
121
122 using Eval = EvaluatorTensorProduct<variant,
123 dim,
124 fe_degree + 1,
125 n_q_points_1d,
126 Number>;
127
128 static void
129 evaluate(const unsigned int n_components,
130 const EvaluationFlags::EvaluationFlags evaluation_flag,
131 const Number * values_dofs_actual,
132 FEEvaluationData<dim, Number, false> & fe_eval);
133
134 static void
135 integrate(const unsigned int n_components,
136 const EvaluationFlags::EvaluationFlags integration_flag,
137 Number * values_dofs_actual,
138 FEEvaluationData<dim, Number, false> & fe_eval,
139 const bool add_into_values_array);
140
141 static Eval
142 create_evaluator_tensor_product(
143 const MatrixFreeFunctions::UnivariateShapeData<Number>
144 *univariate_shape_data)
145 {
146 if (variant == evaluate_evenodd)
147 return Eval(univariate_shape_data->shape_values_eo,
148 univariate_shape_data->shape_gradients_eo,
149 univariate_shape_data->shape_hessians_eo,
150 univariate_shape_data->fe_degree + 1,
151 univariate_shape_data->n_q_points_1d);
152 else
153 return Eval(univariate_shape_data->shape_values,
154 univariate_shape_data->shape_gradients,
155 univariate_shape_data->shape_hessians,
156 univariate_shape_data->fe_degree + 1,
157 univariate_shape_data->n_q_points_1d);
158 }
159 };
160
161
162
167 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
168 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
169 dim,
170 fe_degree,
171 n_q_points_1d,
172 Number>
173 {
174 static void
175 evaluate(const unsigned int n_components,
176 const EvaluationFlags::EvaluationFlags evaluation_flag,
177 const Number * values_dofs_actual,
178 FEEvaluationData<dim, Number, false> & fe_eval);
179
180 static void
181 integrate(const unsigned int n_components,
182 const EvaluationFlags::EvaluationFlags integration_flag,
183 Number * values_dofs_actual,
184 FEEvaluationData<dim, Number, false> & fe_eval,
185 const bool add_into_values_array);
186 };
187
192 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
193 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
194 dim,
195 fe_degree,
196 n_q_points_1d,
197 Number>
198 {
199 template <bool integrate>
200 static void
201 evaluate_or_integrate(
202 const EvaluationFlags::EvaluationFlags evaluation_flag,
203 Number * values_dofs_actual,
204 FEEvaluationData<dim, Number, false> & fe_eval,
205 const bool add_into_values_array = false);
206
207 private:
208 template <typename EvalType>
209 static EvalType
210 create_evaluator_tensor_product(
211 const MatrixFreeFunctions::UnivariateShapeData<Number> &shape_data)
212 {
213 return EvalType(shape_data.shape_values,
214 shape_data.shape_gradients,
215 shape_data.shape_hessians);
216 }
217
218 template <int normal_dir>
219 static void
220 evaluate_tensor_product_per_component(
221 const EvaluationFlags::EvaluationFlags evaluation_flag,
222 Number * values_dofs_actual,
223 FEEvaluationData<dim, Number, false> & fe_eval,
224 const bool add_into_values_array,
225 std::integral_constant<bool, false>);
226
227 template <int normal_dir>
228 static void
229 evaluate_tensor_product_per_component(
230 const EvaluationFlags::EvaluationFlags evaluation_flag,
231 Number * values_dofs_actual,
232 FEEvaluationData<dim, Number, false> & fe_eval,
233 const bool add_into_values_array,
234 std::integral_constant<bool, true>);
235 };
236
237
238
239 template <MatrixFreeFunctions::ElementType type,
240 int dim,
241 int fe_degree,
242 int n_q_points_1d,
243 typename Number>
244 inline void
245 FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::evaluate(
246 const unsigned int n_components,
247 const EvaluationFlags::EvaluationFlags evaluation_flag,
248 const Number * values_dofs_actual,
249 FEEvaluationData<dim, Number, false> & fe_eval)
250 {
251 if (evaluation_flag == EvaluationFlags::nothing)
252 return;
253
254 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
255 univariate_shape_data;
256
257 const auto &shape_data = fe_eval.get_shape_info().data;
258
259 univariate_shape_data.fill(&shape_data.front());
260
261 if (shape_data.size() == dim)
262 for (int i = 1; i < dim; ++i)
263 univariate_shape_data[i] = &shape_data[i];
264
265 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
266 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
267 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
268
269 const unsigned int temp_size =
270 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
271 0 :
272 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
273 Eval::n_rows_of_product :
274 Eval::n_columns_of_product);
275 Number *temp1 = fe_eval.get_scratch_data().begin();
276 Number *temp2;
277 if (temp_size == 0)
278 {
279 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
280 shape_data.front().fe_degree + 1),
281 Utilities::fixed_power<dim>(
282 shape_data.front().n_q_points_1d));
283 }
284 else
285 {
286 temp2 = temp1 + temp_size;
287 }
288
289 const std::size_t n_q_points = temp_size == 0 ?
290 fe_eval.get_shape_info().n_q_points :
291 Eval::n_columns_of_product;
292 const std::size_t dofs_per_comp =
293 (type == MatrixFreeFunctions::truncated_tensor) ?
294 Utilities::pow(shape_data.front().fe_degree + 1, dim) :
295 fe_eval.get_shape_info().dofs_per_component_on_cell;
296 const Number *values_dofs = values_dofs_actual;
297 if (type == MatrixFreeFunctions::truncated_tensor)
298 {
299 const std::size_t n_dofs_per_comp =
300 fe_eval.get_shape_info().dofs_per_component_on_cell;
301 Number *values_dofs_tmp =
302 temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points));
303 const int degree =
304 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
305 for (unsigned int c = 0; c < n_components; ++c)
306 for (int i = 0, count_p = 0, count_q = 0;
307 i < (dim > 2 ? degree + 1 : 1);
308 ++i)
309 {
310 for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
311 {
312 for (int k = 0; k < degree + 1 - j - i;
313 ++k, ++count_p, ++count_q)
314 values_dofs_tmp[c * dofs_per_comp + count_q] =
315 values_dofs_actual[c * n_dofs_per_comp + count_p];
316 for (int k = degree + 1 - j - i; k < degree + 1;
317 ++k, ++count_q)
318 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
319 }
320 for (int j = degree + 1 - i; j < degree + 1; ++j)
321 for (int k = 0; k < degree + 1; ++k, ++count_q)
322 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
323 }
324 values_dofs = values_dofs_tmp;
325 }
326
327 Number *values_quad = fe_eval.begin_values();
328 Number *gradients_quad = fe_eval.begin_gradients();
329 Number *hessians_quad = fe_eval.begin_hessians();
330
331 switch (dim)
332 {
333 case 1:
334 for (unsigned int c = 0; c < n_components; ++c)
335 {
336 if (evaluation_flag & EvaluationFlags::values)
337 eval0.template values<0, true, false>(values_dofs, values_quad);
338 if (evaluation_flag & EvaluationFlags::gradients)
339 eval0.template gradients<0, true, false>(values_dofs,
340 gradients_quad);
341 if (evaluation_flag & EvaluationFlags::hessians)
342 eval0.template hessians<0, true, false>(values_dofs,
343 hessians_quad);
344
345 // advance the next component in 1D array
346 values_dofs += dofs_per_comp;
347 values_quad += n_q_points;
348 gradients_quad += n_q_points;
349 hessians_quad += n_q_points;
350 }
351 break;
352
353 case 2:
354 for (unsigned int c = 0; c < n_components; ++c)
355 {
356 // grad x
357 if (evaluation_flag & EvaluationFlags::gradients)
358 {
359 eval0.template gradients<0, true, false>(values_dofs, temp1);
360 eval1.template values<1, true, false>(temp1, gradients_quad);
361 }
362 if (evaluation_flag & EvaluationFlags::hessians)
363 {
364 // grad xy
365 if (!(evaluation_flag & EvaluationFlags::gradients))
366 eval0.template gradients<0, true, false>(values_dofs,
367 temp1);
368 eval1.template gradients<1, true, false>(temp1,
369 hessians_quad +
370 2 * n_q_points);
371
372 // grad xx
373 eval0.template hessians<0, true, false>(values_dofs, temp1);
374 eval1.template values<1, true, false>(temp1, hessians_quad);
375 }
376
377 // grad y
378 eval0.template values<0, true, false>(values_dofs, temp1);
379 if (evaluation_flag & EvaluationFlags::gradients)
380 eval1.template gradients<1, true, false>(temp1,
381 gradients_quad +
382 n_q_points);
383
384 // grad yy
385 if (evaluation_flag & EvaluationFlags::hessians)
386 eval1.template hessians<1, true, false>(temp1,
387 hessians_quad +
388 n_q_points);
389
390 // val: can use values applied in x
391 if (evaluation_flag & EvaluationFlags::values)
392 eval1.template values<1, true, false>(temp1, values_quad);
393
394 // advance to the next component in 1D array
395 values_dofs += dofs_per_comp;
396 values_quad += n_q_points;
397 gradients_quad += 2 * n_q_points;
398 hessians_quad += 3 * n_q_points;
399 }
400 break;
401
402 case 3:
403 for (unsigned int c = 0; c < n_components; ++c)
404 {
405 if (evaluation_flag & EvaluationFlags::gradients)
406 {
407 // grad x
408 eval0.template gradients<0, true, false>(values_dofs, temp1);
409 eval1.template values<1, true, false>(temp1, temp2);
410 eval2.template values<2, true, false>(temp2, gradients_quad);
411 }
412
413 if (evaluation_flag & EvaluationFlags::hessians)
414 {
415 // grad xz
416 if (!(evaluation_flag & EvaluationFlags::gradients))
417 {
418 eval0.template gradients<0, true, false>(values_dofs,
419 temp1);
420 eval1.template values<1, true, false>(temp1, temp2);
421 }
422 eval2.template gradients<2, true, false>(temp2,
423 hessians_quad +
424 4 * n_q_points);
425
426 // grad xy
427 eval1.template gradients<1, true, false>(temp1, temp2);
428 eval2.template values<2, true, false>(temp2,
429 hessians_quad +
430 3 * n_q_points);
431
432 // grad xx
433 eval0.template hessians<0, true, false>(values_dofs, temp1);
434 eval1.template values<1, true, false>(temp1, temp2);
435 eval2.template values<2, true, false>(temp2, hessians_quad);
436 }
437
438 // grad y
439 eval0.template values<0, true, false>(values_dofs, temp1);
440 if (evaluation_flag & EvaluationFlags::gradients)
441 {
442 eval1.template gradients<1, true, false>(temp1, temp2);
443 eval2.template values<2, true, false>(temp2,
444 gradients_quad +
445 n_q_points);
446 }
447
448 if (evaluation_flag & EvaluationFlags::hessians)
449 {
450 // grad yz
451 if (!(evaluation_flag & EvaluationFlags::gradients))
452 eval1.template gradients<1, true, false>(temp1, temp2);
453 eval2.template gradients<2, true, false>(temp2,
454 hessians_quad +
455 5 * n_q_points);
456
457 // grad yy
458 eval1.template hessians<1, true, false>(temp1, temp2);
459 eval2.template values<2, true, false>(temp2,
460 hessians_quad +
461 n_q_points);
462 }
463
464 // grad z: can use the values applied in x direction stored in
465 // temp1
466 eval1.template values<1, true, false>(temp1, temp2);
467 if (evaluation_flag & EvaluationFlags::gradients)
468 eval2.template gradients<2, true, false>(temp2,
469 gradients_quad +
470 2 * n_q_points);
471
472 // grad zz: can use the values applied in x and y direction stored
473 // in temp2
474 if (evaluation_flag & EvaluationFlags::hessians)
475 eval2.template hessians<2, true, false>(temp2,
476 hessians_quad +
477 2 * n_q_points);
478
479 // val: can use the values applied in x & y direction stored in
480 // temp2
481 if (evaluation_flag & EvaluationFlags::values)
482 eval2.template values<2, true, false>(temp2, values_quad);
483
484 // advance to the next component in 1D array
485 values_dofs += dofs_per_comp;
486 values_quad += n_q_points;
487 gradients_quad += 3 * n_q_points;
488 hessians_quad += 6 * n_q_points;
489 }
490 break;
491
492 default:
493 AssertThrow(false, ExcNotImplemented());
494 }
495
496 // case additional dof for FE_Q_DG0: add values; gradients and second
497 // derivatives evaluate to zero
498 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 &&
499 (evaluation_flag & EvaluationFlags::values))
500 {
501 values_quad -= n_components * n_q_points;
502 values_dofs -= n_components * dofs_per_comp;
503 for (std::size_t c = 0; c < n_components; ++c)
504 for (std::size_t q = 0; q < n_q_points; ++q)
505 values_quad[c * n_q_points + q] +=
506 values_dofs[(c + 1) * dofs_per_comp - 1];
507 }
508 }
509
510
511
512 template <MatrixFreeFunctions::ElementType type,
513 int dim,
514 int fe_degree,
515 int n_q_points_1d,
516 typename Number>
517 inline void
518 FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::integrate(
519 const unsigned int n_components,
520 const EvaluationFlags::EvaluationFlags integration_flag,
521 Number * values_dofs_actual,
522 FEEvaluationData<dim, Number, false> & fe_eval,
523 const bool add_into_values_array)
524 {
525 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
526 univariate_shape_data;
527
528 const auto &shape_data = fe_eval.get_shape_info().data;
529 univariate_shape_data.fill(&shape_data.front());
530
531 if (shape_data.size() == dim)
532 for (int i = 1; i < dim; ++i)
533 univariate_shape_data[i] = &shape_data[i];
534
535 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
536 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
537 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
538
539 const unsigned int temp_size =
540 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
541 0 :
542 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
543 Eval::n_rows_of_product :
544 Eval::n_columns_of_product);
545 Number *temp1 = fe_eval.get_scratch_data().begin();
546 Number *temp2;
547 if (temp_size == 0)
548 {
549 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
550 shape_data.front().fe_degree + 1),
551 Utilities::fixed_power<dim>(
552 shape_data.front().n_q_points_1d));
553 }
554 else
555 {
556 temp2 = temp1 + temp_size;
557 }
558
559 const std::size_t n_q_points = temp_size == 0 ?
560 fe_eval.get_shape_info().n_q_points :
561 Eval::n_columns_of_product;
562 const unsigned int dofs_per_comp =
563 (type == MatrixFreeFunctions::truncated_tensor) ?
564 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
565 fe_eval.get_shape_info().dofs_per_component_on_cell;
566 // expand dof_values to tensor product for truncated tensor products
567 Number *values_dofs =
568 (type == MatrixFreeFunctions::truncated_tensor) ?
569 temp1 + 2 * (std::max<std::size_t>(
570 fe_eval.get_shape_info().dofs_per_component_on_cell,
571 n_q_points)) :
572 values_dofs_actual;
573
574 Number *values_quad = fe_eval.begin_values();
575 Number *gradients_quad = fe_eval.begin_gradients();
576 Number *hessians_quad = fe_eval.begin_hessians();
577
578 switch (dim)
579 {
580 case 1:
581 for (unsigned int c = 0; c < n_components; ++c)
582 {
583 if (integration_flag & EvaluationFlags::values)
584 {
585 if (add_into_values_array == false)
586 eval0.template values<0, false, false>(values_quad,
587 values_dofs);
588 else
589 eval0.template values<0, false, true>(values_quad,
590 values_dofs);
591 }
592 if (integration_flag & EvaluationFlags::gradients)
593 {
594 if (integration_flag & EvaluationFlags::values ||
595 add_into_values_array == true)
596 eval0.template gradients<0, false, true>(gradients_quad,
597 values_dofs);
598 else
599 eval0.template gradients<0, false, false>(gradients_quad,
600 values_dofs);
601 }
602 if ((integration_flag & EvaluationFlags::hessians) != 0u)
603 {
604 if ((integration_flag & EvaluationFlags::values) != 0u ||
605 (integration_flag & EvaluationFlags::gradients) != 0u ||
606 add_into_values_array == true)
607 eval0.template hessians<0, false, true>(hessians_quad,
608 values_dofs);
609 else
610 eval0.template hessians<0, false, false>(hessians_quad,
611 values_dofs);
612 }
613
614 // advance to the next component in 1D array
615 values_dofs += dofs_per_comp;
616 values_quad += n_q_points;
617 gradients_quad += n_q_points;
618 hessians_quad += n_q_points;
619 }
620 break;
621
622 case 2:
623 for (unsigned int c = 0; c < n_components; ++c)
624 {
625 if ((integration_flag & EvaluationFlags::values) &&
626 !(integration_flag & EvaluationFlags::gradients))
627 {
628 eval1.template values<1, false, false>(values_quad, temp1);
629 if (add_into_values_array == false)
630 eval0.template values<0, false, false>(temp1, values_dofs);
631 else
632 eval0.template values<0, false, true>(temp1, values_dofs);
633 }
634 if (integration_flag & EvaluationFlags::gradients)
635 {
636 eval1.template gradients<1, false, false>(gradients_quad +
637 n_q_points,
638 temp1);
639 if (integration_flag & EvaluationFlags::values)
640 eval1.template values<1, false, true>(values_quad, temp1);
641 if (add_into_values_array == false)
642 eval0.template values<0, false, false>(temp1, values_dofs);
643 else
644 eval0.template values<0, false, true>(temp1, values_dofs);
645 eval1.template values<1, false, false>(gradients_quad, temp1);
646 eval0.template gradients<0, false, true>(temp1, values_dofs);
647 }
648 if ((integration_flag & EvaluationFlags::hessians) != 0u)
649 {
650 // grad xx
651 eval1.template values<1, false, false>(hessians_quad, temp1);
652
653 if ((integration_flag & EvaluationFlags::values) != 0u ||
654 (integration_flag & EvaluationFlags::gradients) != 0u ||
655 add_into_values_array == true)
656 eval0.template hessians<0, false, true>(temp1, values_dofs);
657 else
658 eval0.template hessians<0, false, false>(temp1,
659 values_dofs);
660
661 // grad yy
662 eval1.template hessians<1, false, false>(hessians_quad +
663 n_q_points,
664 temp1);
665 eval0.template values<0, false, true>(temp1, values_dofs);
666
667 // grad xy
668 eval1.template gradients<1, false, false>(hessians_quad +
669 2 * n_q_points,
670 temp1);
671 eval0.template gradients<0, false, true>(temp1, values_dofs);
672 }
673
674 // advance to the next component in 1D array
675 values_dofs += dofs_per_comp;
676 values_quad += n_q_points;
677 gradients_quad += 2 * n_q_points;
678 hessians_quad += 3 * n_q_points;
679 }
680 break;
681
682 case 3:
683 for (unsigned int c = 0; c < n_components; ++c)
684 {
685 if ((integration_flag & EvaluationFlags::values) &&
686 !(integration_flag & EvaluationFlags::gradients))
687 {
688 eval2.template values<2, false, false>(values_quad, temp1);
689 eval1.template values<1, false, false>(temp1, temp2);
690 if (add_into_values_array == false)
691 eval0.template values<0, false, false>(temp2, values_dofs);
692 else
693 eval0.template values<0, false, true>(temp2, values_dofs);
694 }
695 if (integration_flag & EvaluationFlags::gradients)
696 {
697 eval2.template gradients<2, false, false>(gradients_quad +
698 2 * n_q_points,
699 temp1);
700 if (integration_flag & EvaluationFlags::values)
701 eval2.template values<2, false, true>(values_quad, temp1);
702 eval1.template values<1, false, false>(temp1, temp2);
703 eval2.template values<2, false, false>(gradients_quad +
704 n_q_points,
705 temp1);
706 eval1.template gradients<1, false, true>(temp1, temp2);
707 if (add_into_values_array == false)
708 eval0.template values<0, false, false>(temp2, values_dofs);
709 else
710 eval0.template values<0, false, true>(temp2, values_dofs);
711 eval2.template values<2, false, false>(gradients_quad, temp1);
712 eval1.template values<1, false, false>(temp1, temp2);
713 eval0.template gradients<0, false, true>(temp2, values_dofs);
714 }
715 if ((integration_flag & EvaluationFlags::hessians) != 0u)
716 {
717 // grad xx
718 eval2.template values<2, false, false>(hessians_quad, temp1);
719 eval1.template values<1, false, false>(temp1, temp2);
720
721 if ((integration_flag & EvaluationFlags::values) != 0u ||
722 (integration_flag & EvaluationFlags::gradients) != 0u ||
723 add_into_values_array == true)
724 eval0.template hessians<0, false, true>(temp2, values_dofs);
725 else
726 eval0.template hessians<0, false, false>(temp2,
727 values_dofs);
728
729 // grad yy
730 eval2.template values<2, false, false>(hessians_quad +
731 n_q_points,
732 temp1);
733 eval1.template hessians<1, false, false>(temp1, temp2);
734 eval0.template values<0, false, true>(temp2, values_dofs);
735
736 // grad zz
737 eval2.template hessians<2, false, false>(hessians_quad +
738 2 * n_q_points,
739 temp1);
740 eval1.template values<1, false, false>(temp1, temp2);
741 eval0.template values<0, false, true>(temp2, values_dofs);
742
743 // grad xy
744 eval2.template values<2, false, false>(hessians_quad +
745 3 * n_q_points,
746 temp1);
747 eval1.template gradients<1, false, false>(temp1, temp2);
748 eval0.template gradients<0, false, true>(temp2, values_dofs);
749
750 // grad xz
751 eval2.template gradients<2, false, false>(hessians_quad +
752 4 * n_q_points,
753 temp1);
754 eval1.template values<1, false, false>(temp1, temp2);
755 eval0.template gradients<0, false, true>(temp2, values_dofs);
756
757 // grad yz
758 eval2.template gradients<2, false, false>(hessians_quad +
759 5 * n_q_points,
760 temp1);
761 eval1.template gradients<1, false, false>(temp1, temp2);
762 eval0.template values<0, false, true>(temp2, values_dofs);
763 }
764
765 // advance to the next component in 1D array
766 values_dofs += dofs_per_comp;
767 values_quad += n_q_points;
768 gradients_quad += 3 * n_q_points;
769 hessians_quad += 6 * n_q_points;
770 }
771 break;
772
773 default:
774 AssertThrow(false, ExcNotImplemented());
775 }
776
777 // case FE_Q_DG0: add values, gradients and second derivatives are zero
778 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
779 {
780 values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
781 values_quad -= n_components * n_q_points;
782 if (integration_flag & EvaluationFlags::values)
783 for (unsigned int c = 0; c < n_components; ++c)
784 {
785 values_dofs[0] = values_quad[0];
786 for (unsigned int q = 1; q < n_q_points; ++q)
787 values_dofs[0] += values_quad[q];
788 values_dofs += dofs_per_comp;
789 values_quad += n_q_points;
790 }
791 else
792 {
793 for (unsigned int c = 0; c < n_components; ++c)
794 values_dofs[c * dofs_per_comp] = Number();
795 values_dofs += n_components * dofs_per_comp;
796 }
797 }
798
799 if (type == MatrixFreeFunctions::truncated_tensor)
800 {
801 const std::size_t n_dofs_per_comp =
802 fe_eval.get_shape_info().dofs_per_component_on_cell;
803 values_dofs -= dofs_per_comp * n_components;
804 const int degree =
805 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
806 for (unsigned int c = 0; c < n_components; ++c)
807 for (int i = 0, count_p = 0, count_q = 0;
808 i < (dim > 2 ? degree + 1 : 1);
809 ++i)
810 {
811 for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
812 {
813 for (int k = 0; k < degree + 1 - j - i;
814 ++k, ++count_p, ++count_q)
815 values_dofs_actual[c * n_dofs_per_comp + count_p] =
816 values_dofs[c * dofs_per_comp + count_q];
817 count_q += j + i;
818 }
819 count_q += i * (degree + 1);
820 }
821 }
822 }
823
824
825
826 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
827 inline void
828 FEEvaluationImpl<
829 MatrixFreeFunctions::tensor_none,
830 dim,
831 fe_degree,
832 n_q_points_1d,
833 Number>::evaluate(const unsigned int n_components,
834 const EvaluationFlags::EvaluationFlags evaluation_flag,
835 const Number * values_dofs_actual,
836 FEEvaluationData<dim, Number, false> & fe_eval)
837 {
838 const std::size_t n_dofs =
839 fe_eval.get_shape_info().dofs_per_component_on_cell;
840 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
841
842 const auto &shape_data = fe_eval.get_shape_info().data;
843
844 using Eval =
845 EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
846
847 if (evaluation_flag & EvaluationFlags::values)
848 {
849 const auto shape_values = shape_data.front().shape_values.data();
850 auto values_quad_ptr = fe_eval.begin_values();
851 auto values_dofs_actual_ptr = values_dofs_actual;
852
853 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
854 for (unsigned int c = 0; c < n_components; ++c)
855 {
856 eval.template values<0, true, false>(values_dofs_actual_ptr,
857 values_quad_ptr);
858
859 values_quad_ptr += n_q_points;
860 values_dofs_actual_ptr += n_dofs;
861 }
862 }
863
864 if (evaluation_flag & EvaluationFlags::gradients)
865 {
866 const auto shape_gradients = shape_data.front().shape_gradients.data();
867 auto gradients_quad_ptr = fe_eval.begin_gradients();
868 auto values_dofs_actual_ptr = values_dofs_actual;
869
870 for (unsigned int c = 0; c < n_components; ++c)
871 {
872 for (unsigned int d = 0; d < dim; ++d)
873 {
874 Eval eval(nullptr,
875 shape_gradients + n_q_points * n_dofs * d,
876 nullptr,
877 n_dofs,
878 n_q_points);
879
880 eval.template gradients<0, true, false>(values_dofs_actual_ptr,
881 gradients_quad_ptr);
882
883 gradients_quad_ptr += n_q_points;
884 }
885 values_dofs_actual_ptr += n_dofs;
886 }
887 }
888
889 if (evaluation_flag & EvaluationFlags::hessians)
890 Assert(false, ExcNotImplemented());
891 }
892
893
894
895 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
896 inline void
897 FEEvaluationImpl<
898 MatrixFreeFunctions::tensor_none,
899 dim,
900 fe_degree,
901 n_q_points_1d,
902 Number>::integrate(const unsigned int n_components,
903 const EvaluationFlags::EvaluationFlags integration_flag,
904 Number * values_dofs_actual,
905 FEEvaluationData<dim, Number, false> &fe_eval,
906 const bool add_into_values_array)
907 {
908 // TODO: implement hessians
909 AssertThrow(!(integration_flag & EvaluationFlags::hessians),
910 ExcNotImplemented());
911
912 const std::size_t n_dofs =
913 fe_eval.get_shape_info().dofs_per_component_on_cell;
914 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
915
916 const auto &shape_data = fe_eval.get_shape_info().data;
917
918 using Eval =
919 EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
920
921 if (integration_flag & EvaluationFlags::values)
922 {
923 const auto shape_values = shape_data.front().shape_values.data();
924 auto values_quad_ptr = fe_eval.begin_values();
925 auto values_dofs_actual_ptr = values_dofs_actual;
926
927 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
928 for (unsigned int c = 0; c < n_components; ++c)
929 {
930 if (add_into_values_array == false)
931 eval.template values<0, false, false>(values_quad_ptr,
932 values_dofs_actual_ptr);
933 else
934 eval.template values<0, false, true>(values_quad_ptr,
935 values_dofs_actual_ptr);
936
937 values_quad_ptr += n_q_points;
938 values_dofs_actual_ptr += n_dofs;
939 }
940 }
941
942 if (integration_flag & EvaluationFlags::gradients)
943 {
944 const auto shape_gradients = shape_data.front().shape_gradients.data();
945 auto gradients_quad_ptr = fe_eval.begin_gradients();
946 auto values_dofs_actual_ptr = values_dofs_actual;
947
948 for (unsigned int c = 0; c < n_components; ++c)
949 {
950 for (unsigned int d = 0; d < dim; ++d)
951 {
952 Eval eval(nullptr,
953 shape_gradients + n_q_points * n_dofs * d,
954 nullptr,
955 n_dofs,
956 n_q_points);
957
958 if ((add_into_values_array == false &&
959 !(integration_flag & EvaluationFlags::values)) &&
960 d == 0)
961 eval.template gradients<0, false, false>(
962 gradients_quad_ptr, values_dofs_actual_ptr);
963 else
964 eval.template gradients<0, false, true>(
965 gradients_quad_ptr, values_dofs_actual_ptr);
966
967 gradients_quad_ptr += n_q_points;
968 }
969 values_dofs_actual_ptr += n_dofs;
970 }
971 }
972 }
973
974
975 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
976 template <bool integrate>
977 inline void
978 FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
979 dim,
980 fe_degree,
981 n_q_points_1d,
982 Number>::
983 evaluate_or_integrate(
984 const EvaluationFlags::EvaluationFlags evaluation_flag,
985 Number * values_dofs_actual,
986 FEEvaluationData<dim, Number, false> & fe_eval,
987 const bool add_into_values_array)
988 {
989 if (evaluation_flag == EvaluationFlags::nothing)
990 return;
991
992 AssertDimension(fe_eval.get_shape_info().data.size(), 2);
993 // First component:
994 evaluate_tensor_product_per_component<0>(
995 evaluation_flag,
996 values_dofs_actual,
997 fe_eval,
998 add_into_values_array,
999 std::integral_constant<bool, integrate>());
1000 // Second component :
1001 evaluate_tensor_product_per_component<1>(
1002 evaluation_flag,
1003 values_dofs_actual,
1004 fe_eval,
1005 add_into_values_array,
1006 std::integral_constant<bool, integrate>());
1007 if (dim == 3)
1008 {
1009 // Third component
1010 evaluate_tensor_product_per_component<2>(
1011 evaluation_flag,
1012 values_dofs_actual,
1013 fe_eval,
1014 add_into_values_array,
1015 std::integral_constant<bool, integrate>());
1016 }
1017 }
1018
1019 // Helper function that applies the 1d evaluation kernels.
1020 // std::integral_constant<bool, false> is the interpolation path, and
1021 // std::integral_constant<bool, true> below is the integration path.
1022 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1023 template <int normal_dir>
1024 inline void
1025 FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
1026 dim,
1027 fe_degree,
1028 n_q_points_1d,
1029 Number>::
1030 evaluate_tensor_product_per_component(
1031 const EvaluationFlags::EvaluationFlags evaluation_flag,
1032 Number * values_dofs_actual,
1033 FEEvaluationData<dim, Number, false> & fe_eval,
1034 const bool add_into_values_array,
1035 std::integral_constant<bool, false>)
1036 {
1037 (void)add_into_values_array;
1038
1039 using EvalNormal =
1040 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
1041 dim,
1042 (fe_degree == -1) ? 1 : fe_degree + 1,
1043 n_q_points_1d,
1044 Number,
1045 normal_dir>;
1046
1047 using EvalTangent =
1048 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
1049 dim,
1050 (fe_degree == -1) ? 1 : fe_degree,
1051 n_q_points_1d,
1052 Number,
1053 normal_dir>;
1054 using Eval0 =
1055 typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1056 using Eval1 =
1057 typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1058 using Eval2 =
1059 typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1060
1061 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
1062 fe_eval.get_shape_info();
1063 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1064 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1065 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1066 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1067 Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1068 ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1069
1070 Number *temp1 = fe_eval.get_scratch_data().begin();
1071 Number *temp2;
1072
1073 temp2 =
1074 temp1 +
1075 std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1076 Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1077
1078 const std::size_t n_q_points = shape_info.n_q_points;
1079 const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1080
1081 // Initial shift depending on component (normal_dir)
1082 Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1083 Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1084 Number *gradients_quad =
1085 fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1086 Number *hessians_quad =
1087 (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1088 fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1089
1090 switch (dim)
1091 {
1092 case 2:
1093 if (evaluation_flag & EvaluationFlags::gradients)
1094 {
1095 eval0.template gradients<0, true, false>(values_dofs, temp1);
1096 eval1.template values<1, true, false>(temp1, gradients_quad);
1097 }
1098 if (evaluation_flag & EvaluationFlags::hessians)
1099 {
1100 // The evaluation/integration here *should* work, however
1101 // the piola transform is not implemented.
1102 AssertThrow(false, ExcNotImplemented());
1103 // grad xy
1104 if (!(evaluation_flag & EvaluationFlags::gradients))
1105 eval0.template gradients<0, true, false>(values_dofs, temp1);
1106 eval1.template gradients<1, true, false>(temp1,
1107 hessians_quad +
1108 2 * n_q_points);
1109
1110 // grad xx
1111 eval0.template hessians<0, true, false>(values_dofs, temp1);
1112 eval1.template values<1, true, false>(temp1, hessians_quad);
1113 }
1114
1115 // grad y
1116 eval0.template values<0, true, false>(values_dofs, temp1);
1117 if (evaluation_flag & EvaluationFlags::gradients)
1118 eval1.template gradients<1, true, false>(temp1,
1119 gradients_quad +
1120 n_q_points);
1121
1122 // grad yy
1123 if (evaluation_flag & EvaluationFlags::hessians)
1124 eval1.template hessians<1, true, false>(temp1,
1125 hessians_quad + n_q_points);
1126
1127 // val: can use values applied in x
1128 if (evaluation_flag & EvaluationFlags::values)
1129 eval1.template values<1, true, false>(temp1, values_quad);
1130 break;
1131 case 3:
1132 if (evaluation_flag & EvaluationFlags::gradients)
1133 {
1134 // grad x
1135 eval0.template gradients<0, true, false>(values_dofs, temp1);
1136 eval1.template values<1, true, false>(temp1, temp2);
1137 eval2.template values<2, true, false>(temp2, gradients_quad);
1138 }
1139
1140 if (evaluation_flag & EvaluationFlags::hessians)
1141 {
1142 // The evaluation/integration here *should* work, however
1143 // the piola transform is not implemented.
1144 AssertThrow(false, ExcNotImplemented());
1145 // grad xz
1146 if (!(evaluation_flag & EvaluationFlags::gradients))
1147 {
1148 eval0.template gradients<0, true, false>(values_dofs, temp1);
1149 eval1.template values<1, true, false>(temp1, temp2);
1150 }
1151 eval2.template gradients<2, true, false>(temp2,
1152 hessians_quad +
1153 4 * n_q_points);
1154
1155 // grad xy
1156 eval1.template gradients<1, true, false>(temp1, temp2);
1157 eval2.template values<2, true, false>(temp2,
1158 hessians_quad +
1159 3 * n_q_points);
1160
1161 // grad xx
1162 eval0.template hessians<0, true, false>(values_dofs, temp1);
1163 eval1.template values<1, true, false>(temp1, temp2);
1164 eval2.template values<2, true, false>(temp2, hessians_quad);
1165 }
1166
1167 // grad y
1168 eval0.template values<0, true, false>(values_dofs, temp1);
1169 if (evaluation_flag & EvaluationFlags::gradients)
1170 {
1171 eval1.template gradients<1, true, false>(temp1, temp2);
1172 eval2.template values<2, true, false>(temp2,
1173 gradients_quad +
1174 n_q_points);
1175 }
1176
1177 if (evaluation_flag & EvaluationFlags::hessians)
1178 {
1179 // grad yz
1180 if (!(evaluation_flag & EvaluationFlags::gradients))
1181 eval1.template gradients<1, true, false>(temp1, temp2);
1182 eval2.template gradients<2, true, false>(temp2,
1183 hessians_quad +
1184 5 * n_q_points);
1185
1186 // grad yy
1187 eval1.template hessians<1, true, false>(temp1, temp2);
1188 eval2.template values<2, true, false>(temp2,
1189 hessians_quad + n_q_points);
1190 }
1191
1192 // grad z: can use the values applied in x direction stored in
1193 // temp1
1194 eval1.template values<1, true, false>(temp1, temp2);
1195 if (evaluation_flag & EvaluationFlags::gradients)
1196 eval2.template gradients<2, true, false>(temp2,
1197 gradients_quad +
1198 2 * n_q_points);
1199
1200 // grad zz: can use the values applied in x and y direction stored
1201 // in temp2
1202 if (evaluation_flag & EvaluationFlags::hessians)
1203 eval2.template hessians<2, true, false>(temp2,
1204 hessians_quad +
1205 2 * n_q_points);
1206
1207 // val: can use the values applied in x & y direction stored in
1208 // temp2
1209 if (evaluation_flag & EvaluationFlags::values)
1210 eval2.template values<2, true, false>(temp2, values_quad);
1211 break;
1212 default:
1213 AssertThrow(false, ExcNotImplemented());
1214 }
1215 }
1216
1217 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1218 template <int normal_dir>
1219 inline void
1220 FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
1221 dim,
1222 fe_degree,
1223 n_q_points_1d,
1224 Number>::
1225 evaluate_tensor_product_per_component(
1226 const EvaluationFlags::EvaluationFlags evaluation_flag,
1227 Number * values_dofs_actual,
1228 FEEvaluationData<dim, Number, false> & fe_eval,
1229 const bool add_into_values_array,
1230 std::integral_constant<bool, true>)
1231 {
1232 using EvalNormal =
1233 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
1234 dim,
1235 (fe_degree == -1) ? 1 : fe_degree + 1,
1236 n_q_points_1d,
1237 Number,
1238 normal_dir>;
1239
1240 using EvalTangent =
1241 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
1242 dim,
1243 (fe_degree == -1) ? 1 : fe_degree,
1244 n_q_points_1d,
1245 Number,
1246 normal_dir>;
1247 using Eval0 =
1248 typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1249 using Eval1 =
1250 typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1251 using Eval2 =
1252 typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1253
1254 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
1255 fe_eval.get_shape_info();
1256 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1257 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1258 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1259 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1260 Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1261 ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1262
1263 Number *temp1 = fe_eval.get_scratch_data().begin();
1264 Number *temp2;
1265
1266 temp2 =
1267 temp1 +
1268 std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1269 Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1270
1271 const std::size_t n_q_points = shape_info.n_q_points;
1272 const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1273
1274 // Initial shift depending on component (normal_dir)
1275 Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1276 Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1277 Number *gradients_quad =
1278 fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1279 Number *hessians_quad =
1280 (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1281 fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1282
1283 // Integrate path
1284 switch (dim)
1285 {
1286 case 2:
1287 if ((evaluation_flag & EvaluationFlags::values) &&
1288 !(evaluation_flag & EvaluationFlags::gradients))
1289 {
1290 eval1.template values<1, false, false>(values_quad, temp1);
1291 if (add_into_values_array == false)
1292 eval0.template values<0, false, false>(temp1, values_dofs);
1293 else
1294 eval0.template values<0, false, true>(temp1, values_dofs);
1295 }
1296 if (evaluation_flag & EvaluationFlags::gradients)
1297 {
1298 eval1.template gradients<1, false, false>(gradients_quad +
1299 n_q_points,
1300 temp1);
1301 if ((evaluation_flag & EvaluationFlags::values))
1302 eval1.template values<1, false, true>(values_quad, temp1);
1303 if (add_into_values_array == false)
1304 eval0.template values<0, false, false>(temp1, values_dofs);
1305 else
1306 eval0.template values<0, false, true>(temp1, values_dofs);
1307 eval1.template values<1, false, false>(gradients_quad, temp1);
1308 eval0.template gradients<0, false, true>(temp1, values_dofs);
1309 }
1310 if (evaluation_flag & EvaluationFlags::hessians)
1311 {
1312 // grad xx
1313 eval1.template values<1, false, false>(hessians_quad, temp1);
1314
1315 if ((evaluation_flag & EvaluationFlags::values) ||
1316 (evaluation_flag & EvaluationFlags::gradients) ||
1317 add_into_values_array == true)
1318 eval0.template hessians<0, false, true>(temp1, values_dofs);
1319 else
1320 eval0.template hessians<0, false, false>(temp1, values_dofs);
1321
1322 // grad yy
1323 eval1.template hessians<1, false, false>(hessians_quad +
1324 n_q_points,
1325 temp1);
1326 eval0.template values<0, false, true>(temp1, values_dofs);
1327
1328 // grad xy
1329 eval1.template gradients<1, false, false>(hessians_quad +
1330 2 * n_q_points,
1331 temp1);
1332 eval0.template gradients<0, false, true>(temp1, values_dofs);
1333 }
1334 break;
1335
1336 case 3:
1337 if ((evaluation_flag & EvaluationFlags::values) &&
1338 !(evaluation_flag & EvaluationFlags::gradients))
1339 {
1340 eval2.template values<2, false, false>(values_quad, temp1);
1341 eval1.template values<1, false, false>(temp1, temp2);
1342 if (add_into_values_array == false)
1343 eval0.template values<0, false, false>(temp2, values_dofs);
1344 else
1345 eval0.template values<0, false, true>(temp2, values_dofs);
1346 }
1347 if (evaluation_flag & EvaluationFlags::gradients)
1348 {
1349 eval2.template gradients<2, false, false>(gradients_quad +
1350 2 * n_q_points,
1351 temp1);
1352 if ((evaluation_flag & EvaluationFlags::values))
1353 eval2.template values<2, false, true>(values_quad, temp1);
1354 eval1.template values<1, false, false>(temp1, temp2);
1355 eval2.template values<2, false, false>(gradients_quad +
1356 n_q_points,
1357 temp1);
1358 eval1.template gradients<1, false, true>(temp1, temp2);
1359 if (add_into_values_array == false)
1360 eval0.template values<0, false, false>(temp2, values_dofs);
1361 else
1362 eval0.template values<0, false, true>(temp2, values_dofs);
1363 eval2.template values<2, false, false>(gradients_quad, temp1);
1364 eval1.template values<1, false, false>(temp1, temp2);
1365 eval0.template gradients<0, false, true>(temp2, values_dofs);
1366 }
1367 if (evaluation_flag & EvaluationFlags::hessians)
1368 {
1369 // grad xx
1370 eval2.template values<2, false, false>(hessians_quad, temp1);
1371 eval1.template values<1, false, false>(temp1, temp2);
1372
1373 if ((evaluation_flag & EvaluationFlags::values) ||
1374 (evaluation_flag & EvaluationFlags::gradients) ||
1375 add_into_values_array == true)
1376 eval0.template hessians<0, false, true>(temp2, values_dofs);
1377 else
1378 eval0.template hessians<0, false, false>(temp2, values_dofs);
1379
1380 // grad yy
1381 eval2.template values<2, false, false>(hessians_quad + n_q_points,
1382 temp1);
1383 eval1.template hessians<1, false, false>(temp1, temp2);
1384 eval0.template values<0, false, true>(temp2, values_dofs);
1385
1386 // grad zz
1387 eval2.template hessians<2, false, false>(hessians_quad +
1388 2 * n_q_points,
1389 temp1);
1390 eval1.template values<1, false, false>(temp1, temp2);
1391 eval0.template values<0, false, true>(temp2, values_dofs);
1392
1393 // grad xy
1394 eval2.template values<2, false, false>(hessians_quad +
1395 3 * n_q_points,
1396 temp1);
1397 eval1.template gradients<1, false, false>(temp1, temp2);
1398 eval0.template gradients<0, false, true>(temp2, values_dofs);
1399
1400 // grad xz
1401 eval2.template gradients<2, false, false>(hessians_quad +
1402 4 * n_q_points,
1403 temp1);
1404 eval1.template values<1, false, false>(temp1, temp2);
1405 eval0.template gradients<0, false, true>(temp2, values_dofs);
1406
1407 // grad yz
1408 eval2.template gradients<2, false, false>(hessians_quad +
1409 5 * n_q_points,
1410 temp1);
1411 eval1.template gradients<1, false, false>(temp1, temp2);
1412 eval0.template values<0, false, true>(temp2, values_dofs);
1413 }
1414
1415 break;
1416 default:
1417 AssertThrow(false, ExcNotImplemented());
1418 }
1419 }
1420
1430 template <EvaluatorVariant variant,
1431 EvaluatorQuantity quantity,
1432 int dim,
1433 int basis_size_1,
1434 int basis_size_2,
1435 typename Number,
1436 typename Number2>
1437 struct FEEvaluationImplBasisChange
1438 {
1439 static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
1440 "The second dimension must not be smaller than the first");
1441
1464#ifndef DEBUG
1465 DEAL_II_ALWAYS_INLINE
1466#endif
1467 static void
1468 do_forward(
1469 const unsigned int n_components,
1470 const AlignedVector<Number2> &transformation_matrix,
1471 const Number * values_in,
1472 Number * values_out,
1473 const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1474 const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1475 {
1476 Assert(
1477 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1478 ExcMessage("The second dimension must not be smaller than the first"));
1479
1480 Assert(quantity == EvaluatorQuantity::value, ExcInternalError());
1481
1482 // we do recursion until dim==1 or dim==2 and we have
1483 // basis_size_1==basis_size_2. The latter optimization increases
1484 // optimization possibilities for the compiler but does only work for
1485 // aliased pointers if the sizes are equal.
1486 constexpr int next_dim =
1487 (dim > 2 ||
1488 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1489 dim - 1 :
1490 dim;
1491
1492 EvaluatorTensorProduct<variant,
1493 dim,
1494 basis_size_1,
1495 (basis_size_1 == 0 ? 0 : basis_size_2),
1496 Number,
1497 Number2>
1498 eval_val(transformation_matrix,
1499 AlignedVector<Number2>(),
1500 AlignedVector<Number2>(),
1501 basis_size_1_variable,
1502 basis_size_2_variable);
1503 const unsigned int np_1 =
1504 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1505 const unsigned int np_2 =
1506 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1507 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1508 ExcMessage("Cannot transform with 0-point basis"));
1509 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1510 ExcMessage("Cannot transform with 0-point basis"));
1511
1512 // run loop backwards to ensure correctness if values_in aliases with
1513 // values_out in case with basis_size_1 < basis_size_2
1514 values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
1515 values_out =
1516 values_out + n_components * Utilities::fixed_power<dim>(np_2);
1517 for (unsigned int c = n_components; c != 0; --c)
1518 {
1519 values_in -= Utilities::fixed_power<dim>(np_1);
1520 values_out -= Utilities::fixed_power<dim>(np_2);
1521 if (next_dim < dim)
1522 for (unsigned int q = np_1; q != 0; --q)
1523 FEEvaluationImplBasisChange<
1524 variant,
1525 quantity,
1526 next_dim,
1527 basis_size_1,
1528 basis_size_2,
1529 Number,
1530 Number2>::do_forward(1,
1531 transformation_matrix,
1532 values_in +
1533 (q - 1) *
1534 Utilities::fixed_power<next_dim>(np_1),
1535 values_out +
1536 (q - 1) *
1537 Utilities::fixed_power<next_dim>(np_2),
1538 basis_size_1_variable,
1539 basis_size_2_variable);
1540
1541 // the recursion stops if dim==1 or if dim==2 and
1542 // basis_size_1==basis_size_2 (the latter is used because the
1543 // compiler generates nicer code)
1544 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1545 {
1546 eval_val.template values<0, true, false>(values_in, values_out);
1547 eval_val.template values<1, true, false>(values_out, values_out);
1548 }
1549 else if (dim == 1)
1550 eval_val.template values<dim - 1, true, false>(values_in,
1551 values_out);
1552 else
1553 eval_val.template values<dim - 1, true, false>(values_out,
1554 values_out);
1555 }
1556 }
1557
1588#ifndef DEBUG
1589 DEAL_II_ALWAYS_INLINE
1590#endif
1591 static void
1592 do_backward(
1593 const unsigned int n_components,
1594 const AlignedVector<Number2> &transformation_matrix,
1595 const bool add_into_result,
1596 Number * values_in,
1597 Number * values_out,
1598 const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1599 const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1600 {
1601 Assert(
1602 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1603 ExcMessage("The second dimension must not be smaller than the first"));
1604 Assert(add_into_result == false || values_in != values_out,
1605 ExcMessage(
1606 "Input and output cannot alias with each other when "
1607 "adding the result of the basis change to existing data"));
1608
1609 Assert(quantity == EvaluatorQuantity::value ||
1610 quantity == EvaluatorQuantity::hessian,
1611 ExcInternalError());
1612
1613 constexpr int next_dim =
1614 (dim > 2 ||
1615 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1616 dim - 1 :
1617 dim;
1618 EvaluatorTensorProduct<variant,
1619 dim,
1620 basis_size_1,
1621 (basis_size_1 == 0 ? 0 : basis_size_2),
1622 Number,
1623 Number2>
1624 eval_val(transformation_matrix,
1625 transformation_matrix,
1626 transformation_matrix,
1627 basis_size_1_variable,
1628 basis_size_2_variable);
1629 const unsigned int np_1 =
1630 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1631 const unsigned int np_2 =
1632 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1633 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1634 ExcMessage("Cannot transform with 0-point basis"));
1635 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1636 ExcMessage("Cannot transform with 0-point basis"));
1637
1638 for (unsigned int c = 0; c < n_components; ++c)
1639 {
1640 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1641 {
1642 if (quantity == EvaluatorQuantity::value)
1643 eval_val.template values<1, false, false>(values_in, values_in);
1644 else
1645 eval_val.template hessians<1, false, false>(values_in,
1646 values_in);
1647
1648 if (add_into_result)
1649 {
1650 if (quantity == EvaluatorQuantity::value)
1651 eval_val.template values<0, false, true>(values_in,
1652 values_out);
1653 else
1654 eval_val.template hessians<0, false, true>(values_in,
1655 values_out);
1656 }
1657 else
1658 {
1659 if (quantity == EvaluatorQuantity::value)
1660 eval_val.template values<0, false, false>(values_in,
1661 values_out);
1662 else
1663 eval_val.template hessians<0, false, false>(values_in,
1664 values_out);
1665 }
1666 }
1667 else
1668 {
1669 if (dim == 1 && add_into_result)
1670 {
1671 if (quantity == EvaluatorQuantity::value)
1672 eval_val.template values<0, false, true>(values_in,
1673 values_out);
1674 else
1675 eval_val.template hessians<0, false, true>(values_in,
1676 values_out);
1677 }
1678 else if (dim == 1)
1679 {
1680 if (quantity == EvaluatorQuantity::value)
1681 eval_val.template values<0, false, false>(values_in,
1682 values_out);
1683 else
1684 eval_val.template hessians<0, false, false>(values_in,
1685 values_out);
1686 }
1687 else
1688 {
1689 if (quantity == EvaluatorQuantity::value)
1690 eval_val.template values<dim - 1, false, false>(values_in,
1691 values_in);
1692 else
1693 eval_val.template hessians<dim - 1, false, false>(
1694 values_in, values_in);
1695 }
1696 }
1697 if (next_dim < dim)
1698 for (unsigned int q = 0; q < np_1; ++q)
1699 FEEvaluationImplBasisChange<variant,
1700 quantity,
1701 next_dim,
1702 basis_size_1,
1703 basis_size_2,
1704 Number,
1705 Number2>::
1706 do_backward(1,
1707 transformation_matrix,
1708 add_into_result,
1709 values_in +
1710 q * Utilities::fixed_power<next_dim>(np_2),
1711 values_out +
1712 q * Utilities::fixed_power<next_dim>(np_1),
1713 basis_size_1_variable,
1714 basis_size_2_variable);
1715
1716 values_in += Utilities::fixed_power<dim>(np_2);
1717 values_out += Utilities::fixed_power<dim>(np_1);
1718 }
1719 }
1720
1741 static void
1742 do_mass(const unsigned int n_components,
1743 const AlignedVector<Number2> &transformation_matrix,
1744 const AlignedVector<Number> & coefficients,
1745 const Number * values_in,
1746 Number * scratch_data,
1747 Number * values_out)
1748 {
1749 constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1750 Number * my_scratch =
1751 basis_size_1 != basis_size_2 ? scratch_data : values_out;
1752
1753 const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1754 Assert(coefficients.size() == size_per_component ||
1755 coefficients.size() == n_components * size_per_component,
1756 ExcDimensionMismatch(coefficients.size(), size_per_component));
1757 const unsigned int stride =
1758 coefficients.size() == size_per_component ? 0 : 1;
1759
1760 for (unsigned int q = basis_size_1; q != 0; --q)
1761 FEEvaluationImplBasisChange<
1762 variant,
1763 EvaluatorQuantity::value,
1764 next_dim,
1765 basis_size_1,
1766 basis_size_2,
1767 Number,
1768 Number2>::do_forward(n_components,
1769 transformation_matrix,
1770 values_in +
1771 (q - 1) *
1772 Utilities::pow(basis_size_1, dim - 1),
1773 my_scratch +
1774 (q - 1) *
1775 Utilities::pow(basis_size_2, dim - 1));
1776 EvaluatorTensorProduct<variant,
1777 dim,
1778 basis_size_1,
1779 basis_size_2,
1780 Number,
1781 Number2>
1782 eval_val(transformation_matrix);
1783 const unsigned int n_inner_blocks =
1784 (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1785 const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1786 for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1787 for (unsigned int c = 0; c < n_components; ++c)
1788 {
1789 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1790 eval_val.template values_one_line<dim - 1, true, false>(
1791 my_scratch + i, my_scratch + i);
1792 for (unsigned int q = 0; q < basis_size_2; ++q)
1793 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1794 my_scratch[i + q * n_blocks + c * size_per_component] *=
1795 coefficients[i + q * n_blocks +
1796 c * stride * size_per_component];
1797 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1798 eval_val.template values_one_line<dim - 1, false, false>(
1799 my_scratch + i, my_scratch + i);
1800 }
1801 for (unsigned int q = 0; q < basis_size_1; ++q)
1802 FEEvaluationImplBasisChange<
1803 variant,
1804 EvaluatorQuantity::value,
1805 next_dim,
1806 basis_size_1,
1807 basis_size_2,
1808 Number,
1809 Number2>::do_backward(n_components,
1810 transformation_matrix,
1811 false,
1812 my_scratch +
1813 q * Utilities::pow(basis_size_2, dim - 1),
1814 values_out +
1815 q * Utilities::pow(basis_size_1, dim - 1));
1816 }
1817 };
1818
1819
1820
1833 template <int dim, int fe_degree, typename Number>
1834 struct FEEvaluationImplCollocation
1835 {
1836 static void
1837 evaluate(const unsigned int n_components,
1838 const EvaluationFlags::EvaluationFlags evaluation_flag,
1839 const Number * values_dofs,
1840 FEEvaluationData<dim, Number, false> & fe_eval);
1841
1842 static void
1843 do_evaluate(const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
1844 const EvaluationFlags::EvaluationFlags evaluation_flag,
1845 const Number * values_dofs,
1846 Number * gradients_quad,
1847 Number * hessians_quad);
1848
1849 static void
1850 integrate(const unsigned int n_components,
1851 const EvaluationFlags::EvaluationFlags integration_flag,
1852 Number * values_dofs,
1853 FEEvaluationData<dim, Number, false> & fe_eval,
1854 const bool add_into_values_array);
1855
1856 static void
1857 do_integrate(const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
1858 const EvaluationFlags::EvaluationFlags integration_flag,
1859 Number * values_dofs,
1860 Number * gradients_quad,
1861 const Number * hessians_quad,
1862 const bool add_into_values_array);
1863 };
1864
1865
1866
1867 template <int dim, int fe_degree, typename Number>
1868 inline void
1869 FEEvaluationImplCollocation<dim, fe_degree, Number>::evaluate(
1870 const unsigned int n_components,
1871 const EvaluationFlags::EvaluationFlags evaluation_flag,
1872 const Number * values_dofs,
1873 FEEvaluationData<dim, Number, false> & fe_eval)
1874 {
1875 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1876
1877 for (unsigned int c = 0; c < n_components; ++c)
1878 {
1879 if ((evaluation_flag & EvaluationFlags::values) != 0u)
1880 for (unsigned int i = 0; i < n_points; ++i)
1881 fe_eval.begin_values()[n_points * c + i] =
1882 values_dofs[n_points * c + i];
1883
1884 do_evaluate(fe_eval.get_shape_info().data.front(),
1885 evaluation_flag,
1886 values_dofs + c * n_points,
1887 fe_eval.begin_gradients() + c * dim * n_points,
1888 fe_eval.begin_hessians() +
1889 c * dim * (dim + 1) / 2 * n_points);
1890 }
1891 }
1892
1893
1894
1895 template <int dim, int fe_degree, typename Number>
1896 inline void
1897 FEEvaluationImplCollocation<dim, fe_degree, Number>::do_evaluate(
1898 const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
1899 const EvaluationFlags::EvaluationFlags evaluation_flag,
1900 const Number * values_dofs,
1901 Number * gradients_quad,
1902 Number * hessians_quad)
1903 {
1904 AssertDimension(shape.shape_gradients_collocation_eo.size(),
1905 (fe_degree + 2) / 2 * (fe_degree + 1));
1906 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1907
1908 EvaluatorTensorProduct<evaluate_evenodd,
1909 dim,
1910 fe_degree + 1,
1911 fe_degree + 1,
1912 Number>
1913 eval(AlignedVector<Number>(),
1914 shape.shape_gradients_collocation_eo,
1915 shape.shape_hessians_collocation_eo);
1916 if ((evaluation_flag &
1917 (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u)
1918 {
1919 eval.template gradients<0, true, false>(values_dofs, gradients_quad);
1920 if (dim > 1)
1921 eval.template gradients<1, true, false>(values_dofs,
1922 gradients_quad + n_points);
1923 if (dim > 2)
1924 eval.template gradients<2, true, false>(values_dofs,
1925 gradients_quad +
1926 2 * n_points);
1927 }
1928 if (evaluation_flag & EvaluationFlags::hessians)
1929 {
1930 eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1931 if (dim > 1)
1932 {
1933 eval.template gradients<1, true, false>(gradients_quad,
1934 hessians_quad +
1935 dim * n_points);
1936 eval.template hessians<1, true, false>(values_dofs,
1937 hessians_quad + n_points);
1938 }
1939 if (dim > 2)
1940 {
1941 eval.template gradients<2, true, false>(gradients_quad,
1942 hessians_quad +
1943 4 * n_points);
1944 eval.template gradients<2, true, false>(gradients_quad + n_points,
1945 hessians_quad +
1946 5 * n_points);
1947 eval.template hessians<2, true, false>(values_dofs,
1948 hessians_quad +
1949 2 * n_points);
1950 }
1951 }
1952 }
1953
1954
1955
1956 template <int dim, int fe_degree, typename Number>
1957 inline void
1958 FEEvaluationImplCollocation<dim, fe_degree, Number>::integrate(
1959 const unsigned int n_components,
1960 const EvaluationFlags::EvaluationFlags integration_flag,
1961 Number * values_dofs,
1962 FEEvaluationData<dim, Number, false> & fe_eval,
1963 const bool add_into_values_array)
1964 {
1965 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1966
1967 for (unsigned int c = 0; c < n_components; ++c)
1968 {
1969 if ((integration_flag & EvaluationFlags::values) != 0u)
1970 {
1971 if (add_into_values_array)
1972 for (unsigned int i = 0; i < n_points; ++i)
1973 values_dofs[n_points * c + i] +=
1974 fe_eval.begin_values()[n_points * c + i];
1975 else
1976 for (unsigned int i = 0; i < n_points; ++i)
1977 values_dofs[n_points * c + i] =
1978 fe_eval.begin_values()[n_points * c + i];
1979 }
1980
1981 do_integrate(fe_eval.get_shape_info().data.front(),
1982 integration_flag,
1983 values_dofs + c * n_points,
1984 fe_eval.begin_gradients() + c * dim * n_points,
1985 fe_eval.begin_hessians() +
1986 c * dim * (dim + 1) / 2 * n_points,
1987 add_into_values_array ||
1988 ((integration_flag & EvaluationFlags::values) != 0u));
1989 }
1990 }
1991
1992
1993
1994 template <int dim, int fe_degree, typename Number>
1995 inline void
1996 FEEvaluationImplCollocation<dim, fe_degree, Number>::do_integrate(
1997 const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
1998 const EvaluationFlags::EvaluationFlags integration_flag,
1999 Number * values_dofs,
2000 Number * gradients_quad,
2001 const Number * hessians_quad,
2002 const bool add_into_values_array)
2003 {
2004 AssertDimension(shape.shape_gradients_collocation_eo.size(),
2005 (fe_degree + 2) / 2 * (fe_degree + 1));
2006
2007 EvaluatorTensorProduct<evaluate_evenodd,
2008 dim,
2009 fe_degree + 1,
2010 fe_degree + 1,
2011 Number>
2012 eval(AlignedVector<Number>(),
2013 shape.shape_gradients_collocation_eo,
2014 shape.shape_hessians_collocation_eo);
2015 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
2016
2017 if ((integration_flag & EvaluationFlags::hessians) != 0u)
2018 {
2019 // diagonal
2020 // grad xx
2021 if (add_into_values_array == true)
2022 eval.template hessians<0, false, true>(hessians_quad, values_dofs);
2023 else
2024 eval.template hessians<0, false, false>(hessians_quad, values_dofs);
2025 // grad yy
2026 if (dim > 1)
2027 eval.template hessians<1, false, true>(hessians_quad + n_points,
2028 values_dofs);
2029 // grad zz
2030 if (dim > 2)
2031 eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
2032 values_dofs);
2033 // off-diagonal
2034 if (dim == 2)
2035 {
2036 // grad xy, queue into gradient
2037 if (integration_flag & EvaluationFlags::gradients)
2038 eval.template gradients<1, false, true>(hessians_quad +
2039 2 * n_points,
2040 gradients_quad);
2041 else
2042 eval.template gradients<1, false, false>(hessians_quad +
2043 2 * n_points,
2044 gradients_quad);
2045 }
2046 if (dim == 3)
2047 {
2048 // grad xy, queue into gradient
2049 if (integration_flag & EvaluationFlags::gradients)
2050 eval.template gradients<1, false, true>(hessians_quad +
2051 3 * n_points,
2052 gradients_quad);
2053 else
2054 eval.template gradients<1, false, false>(hessians_quad +
2055 3 * n_points,
2056 gradients_quad);
2057
2058 // grad xz
2059 eval.template gradients<2, false, true>(hessians_quad +
2060 4 * n_points,
2061 gradients_quad);
2062
2063 // grad yz
2064 if (integration_flag & EvaluationFlags::gradients)
2065 eval.template gradients<2, false, true>(
2066 hessians_quad + 5 * n_points, gradients_quad + n_points);
2067 else
2068 eval.template gradients<2, false, false>(
2069 hessians_quad + 5 * n_points, gradients_quad + n_points);
2070 }
2071
2072 // if we did not integrate gradients, set the last slot to zero
2073 // which was not touched before, in order to avoid the if
2074 // statement in the gradients loop below
2075 if ((integration_flag & EvaluationFlags::gradients) == 0u)
2076 for (unsigned int q = 0; q < n_points; ++q)
2077 gradients_quad[(dim - 1) * n_points + q] = Number();
2078 }
2079
2080 if ((integration_flag &
2081 (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u)
2082 {
2083 if (add_into_values_array ||
2084 (integration_flag & EvaluationFlags::hessians) != 0u)
2085 eval.template gradients<0, false, true>(gradients_quad, values_dofs);
2086 else
2087 eval.template gradients<0, false, false>(gradients_quad, values_dofs);
2088 if (dim > 1)
2089 eval.template gradients<1, false, true>(gradients_quad + n_points,
2090 values_dofs);
2091 if (dim > 2)
2092 eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
2093 values_dofs);
2094 }
2095 }
2096
2097
2098
2109 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2110 struct FEEvaluationImplTransformToCollocation
2111 {
2112 static void
2113 evaluate(const unsigned int n_components,
2114 const EvaluationFlags::EvaluationFlags evaluation_flag,
2115 const Number * values_dofs,
2116 FEEvaluationData<dim, Number, false> & fe_eval);
2117
2118 static void
2119 integrate(const unsigned int n_components,
2120 const EvaluationFlags::EvaluationFlags evaluation_flag,
2121 Number * values_dofs,
2122 FEEvaluationData<dim, Number, false> & fe_eval,
2123 const bool add_into_values_array);
2124 };
2125
2126
2127
2128 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2129 inline void
2130 FEEvaluationImplTransformToCollocation<
2131 dim,
2132 fe_degree,
2133 n_q_points_1d,
2134 Number>::evaluate(const unsigned int n_components,
2135 const EvaluationFlags::EvaluationFlags evaluation_flag,
2136 const Number * values_dofs,
2137 FEEvaluationData<dim, Number, false> & fe_eval)
2138 {
2139 const auto &shape_data = fe_eval.get_shape_info().data.front();
2140
2141 Assert(n_q_points_1d > fe_degree,
2142 ExcMessage("You lose information when going to a collocation space "
2143 "of lower degree, so the evaluation results would be "
2144 "wrong. Thus, this class does not permit the desired "
2145 "operation."));
2146 constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
2147 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2148
2149 for (unsigned int c = 0; c < n_components; ++c)
2150 {
2151 FEEvaluationImplBasisChange<
2152 evaluate_evenodd,
2153 EvaluatorQuantity::value,
2154 dim,
2155 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2156 n_q_points_1d,
2157 Number,
2158 Number>::do_forward(1,
2159 shape_data.shape_values_eo,
2160 values_dofs + c * n_dofs,
2161 fe_eval.begin_values() + c * n_q_points);
2162
2163 // apply derivatives in the collocation space
2164 if (evaluation_flag &
2165 (EvaluationFlags::gradients | EvaluationFlags::hessians))
2166 FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
2167 do_evaluate(shape_data,
2168 evaluation_flag & (EvaluationFlags::gradients |
2169 EvaluationFlags::hessians),
2170 fe_eval.begin_values() + c * n_q_points,
2171 fe_eval.begin_gradients() + c * dim * n_q_points,
2172 fe_eval.begin_hessians() +
2173 c * dim * (dim + 1) / 2 * n_q_points);
2174 }
2175 }
2176
2177
2178
2179 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2180 inline void
2181 FEEvaluationImplTransformToCollocation<
2182 dim,
2183 fe_degree,
2184 n_q_points_1d,
2185 Number>::integrate(const unsigned int n_components,
2186 const EvaluationFlags::EvaluationFlags integration_flag,
2187 Number * values_dofs,
2188 FEEvaluationData<dim, Number, false> & fe_eval,
2189 const bool add_into_values_array)
2190 {
2191 const auto &shape_data = fe_eval.get_shape_info().data.front();
2192
2193 Assert(n_q_points_1d > fe_degree,
2194 ExcMessage("You lose information when going to a collocation space "
2195 "of lower degree, so the evaluation results would be "
2196 "wrong. Thus, this class does not permit the desired "
2197 "operation."));
2198 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2199
2200 for (unsigned int c = 0; c < n_components; ++c)
2201 {
2202 // apply derivatives in collocation space
2203 if (integration_flag &
2204 (EvaluationFlags::gradients | EvaluationFlags::hessians))
2205 FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
2206 do_integrate(shape_data,
2207 integration_flag & (EvaluationFlags::gradients |
2208 EvaluationFlags::hessians),
2209 fe_eval.begin_values() + c * n_q_points,
2210 fe_eval.begin_gradients() + c * dim * n_q_points,
2211 fe_eval.begin_hessians() +
2212 c * dim * (dim + 1) / 2 * n_q_points,
2213 /*add_into_values_array=*/
2214 integration_flag & EvaluationFlags::values);
2215
2216 // transform back to the original space
2217 FEEvaluationImplBasisChange<
2218 evaluate_evenodd,
2219 EvaluatorQuantity::value,
2220 dim,
2221 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2222 n_q_points_1d,
2223 Number,
2224 Number>::do_backward(1,
2225 shape_data.shape_values_eo,
2226 add_into_values_array,
2227 fe_eval.begin_values() + c * n_q_points,
2228 values_dofs +
2229 c * Utilities::pow(fe_degree + 1, dim));
2230 }
2231 }
2232
2233
2234
2242 constexpr bool
2243 use_collocation_evaluation(const unsigned int fe_degree,
2244 const unsigned int n_q_points_1d)
2245 {
2246 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2247 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2248 }
2249
2250
2266 template <int dim, typename Number>
2267 struct FEEvaluationImplEvaluateSelector
2268 {
2269 template <int fe_degree, int n_q_points_1d>
2270 static bool
2271 run(const unsigned int n_components,
2272 const EvaluationFlags::EvaluationFlags evaluation_flag,
2273 const Number * values_dofs,
2274 FEEvaluationData<dim, Number, false> & fe_eval)
2275 {
2276 const auto element_type = fe_eval.get_shape_info().element_type;
2277 using ElementType = MatrixFreeFunctions::ElementType;
2278
2279 Assert(fe_eval.get_shape_info().data.size() == 1 ||
2280 (fe_eval.get_shape_info().data.size() == dim &&
2281 element_type == ElementType::tensor_general) ||
2282 element_type == ElementType::tensor_raviart_thomas,
2283 ExcNotImplemented());
2284
2285 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2286 element_type == ElementType::tensor_symmetric_collocation)
2287 {
2288 FEEvaluationImplCollocation<dim, fe_degree, Number>::evaluate(
2289 n_components, evaluation_flag, values_dofs, fe_eval);
2290 }
2291 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2292 // shape_info.h for more details
2293 else if (fe_degree >= 0 &&
2294 use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2295 element_type <= ElementType::tensor_symmetric)
2296 {
2297 FEEvaluationImplTransformToCollocation<
2298 dim,
2299 fe_degree,
2300 n_q_points_1d,
2301 Number>::evaluate(n_components,
2302 evaluation_flag,
2303 values_dofs,
2304 fe_eval);
2305 }
2306 else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
2307 {
2308 FEEvaluationImpl<ElementType::tensor_symmetric,
2309 dim,
2310 fe_degree,
2311 n_q_points_1d,
2312 Number>::evaluate(n_components,
2313 evaluation_flag,
2314 values_dofs,
2315 fe_eval);
2316 }
2317 else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2318 {
2319 FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
2320 dim,
2321 fe_degree,
2322 n_q_points_1d,
2323 Number>::evaluate(n_components,
2324 evaluation_flag,
2325 values_dofs,
2326 fe_eval);
2327 }
2328 else if (element_type == ElementType::truncated_tensor)
2329 {
2330 FEEvaluationImpl<ElementType::truncated_tensor,
2331 dim,
2332 fe_degree,
2333 n_q_points_1d,
2334 Number>::evaluate(n_components,
2335 evaluation_flag,
2336 values_dofs,
2337 fe_eval);
2338 }
2339 else if (element_type == ElementType::tensor_none)
2340 {
2341 FEEvaluationImpl<ElementType::tensor_none,
2342 dim,
2343 fe_degree,
2344 n_q_points_1d,
2345 Number>::evaluate(n_components,
2346 evaluation_flag,
2347 values_dofs,
2348 fe_eval);
2349 }
2350 else if (element_type == ElementType::tensor_raviart_thomas)
2351 {
2352 FEEvaluationImpl<
2353 ElementType::tensor_raviart_thomas,
2354 dim,
2355 (fe_degree == -1) ? 1 : fe_degree,
2356 (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2357 Number>::template evaluate_or_integrate<false>(evaluation_flag,
2358 const_cast<Number *>(
2359 values_dofs),
2360 fe_eval);
2361 }
2362 else
2363 {
2364 FEEvaluationImpl<ElementType::tensor_general,
2365 dim,
2366 fe_degree,
2367 n_q_points_1d,
2368 Number>::evaluate(n_components,
2369 evaluation_flag,
2370 values_dofs,
2371 fe_eval);
2372 }
2373
2374 return false;
2375 }
2376 };
2377
2378
2379
2395 template <int dim, typename Number>
2396 struct FEEvaluationImplIntegrateSelector
2397 {
2398 template <int fe_degree, int n_q_points_1d>
2399 static bool
2400 run(const unsigned int n_components,
2401 const EvaluationFlags::EvaluationFlags integration_flag,
2402 Number * values_dofs,
2403 FEEvaluationData<dim, Number, false> & fe_eval,
2404 const bool sum_into_values_array)
2405 {
2406 const auto element_type = fe_eval.get_shape_info().element_type;
2407 using ElementType = MatrixFreeFunctions::ElementType;
2408
2409 Assert(fe_eval.get_shape_info().data.size() == 1 ||
2410 (fe_eval.get_shape_info().data.size() == dim &&
2411 element_type == ElementType::tensor_general) ||
2412 element_type == ElementType::tensor_raviart_thomas,
2413 ExcNotImplemented());
2414
2415 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2416 element_type == ElementType::tensor_symmetric_collocation)
2417 {
2418 FEEvaluationImplCollocation<dim, fe_degree, Number>::integrate(
2419 n_components,
2420 integration_flag,
2421 values_dofs,
2422 fe_eval,
2423 sum_into_values_array);
2424 }
2425 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2426 // shape_info.h for more details
2427 else if (fe_degree >= 0 &&
2428 use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2429 element_type <= ElementType::tensor_symmetric)
2430 {
2431 FEEvaluationImplTransformToCollocation<
2432 dim,
2433 fe_degree,
2434 n_q_points_1d,
2435 Number>::integrate(n_components,
2436 integration_flag,
2437 values_dofs,
2438 fe_eval,
2439 sum_into_values_array);
2440 }
2441 else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
2442 {
2443 FEEvaluationImpl<ElementType::tensor_symmetric,
2444 dim,
2445 fe_degree,
2446 n_q_points_1d,
2447 Number>::integrate(n_components,
2448 integration_flag,
2449 values_dofs,
2450 fe_eval,
2451 sum_into_values_array);
2452 }
2453 else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2454 {
2455 FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
2456 dim,
2457 fe_degree,
2458 n_q_points_1d,
2459 Number>::integrate(n_components,
2460 integration_flag,
2461 values_dofs,
2462 fe_eval,
2463 sum_into_values_array);
2464 }
2465 else if (element_type == ElementType::truncated_tensor)
2466 {
2467 FEEvaluationImpl<ElementType::truncated_tensor,
2468 dim,
2469 fe_degree,
2470 n_q_points_1d,
2471 Number>::integrate(n_components,
2472 integration_flag,
2473 values_dofs,
2474 fe_eval,
2475 sum_into_values_array);
2476 }
2477 else if (element_type == ElementType::tensor_none)
2478 {
2479 FEEvaluationImpl<ElementType::tensor_none,
2480 dim,
2481 fe_degree,
2482 n_q_points_1d,
2483 Number>::integrate(n_components,
2484 integration_flag,
2485 values_dofs,
2486 fe_eval,
2487 sum_into_values_array);
2488 }
2489 else if (element_type == ElementType::tensor_raviart_thomas)
2490 {
2491 FEEvaluationImpl<ElementType::tensor_raviart_thomas,
2492 dim,
2493 (fe_degree == -1) ? 1 : fe_degree,
2494 (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2495 Number>::
2496 template evaluate_or_integrate<true>(integration_flag,
2497 const_cast<Number *>(
2498 values_dofs),
2499 fe_eval,
2500 sum_into_values_array);
2501 }
2502 else
2503 {
2504 FEEvaluationImpl<ElementType::tensor_general,
2505 dim,
2506 fe_degree,
2507 n_q_points_1d,
2508 Number>::integrate(n_components,
2509 integration_flag,
2510 values_dofs,
2511 fe_eval,
2512 sum_into_values_array);
2513 }
2514
2515 return false;
2516 }
2517 };
2518
2519
2520
2521 template <bool symmetric_evaluate,
2522 int dim,
2523 int fe_degree,
2524 int n_q_points_1d,
2525 typename Number>
2526 struct FEFaceEvaluationImpl
2527 {
2528 // We enable a transformation to collocation for derivatives if it gives
2529 // correct results (first two conditions), if it is the most efficient
2530 // choice in terms of operation counts (third condition) and if we were
2531 // able to initialize the fields in shape_info.templates.h from the
2532 // polynomials (fourth condition).
2533 using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
2534 evaluate_general,
2535 dim - 1,
2536 fe_degree + 1,
2537 n_q_points_1d,
2538 Number>;
2539
2540 static Eval
2541 create_evaluator_tensor_product(
2542 const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
2543 const unsigned int subface_index,
2544 const unsigned int direction)
2545 {
2546 if (symmetric_evaluate)
2547 return Eval(data.shape_values_eo,
2548 data.shape_gradients_eo,
2549 data.shape_hessians_eo,
2550 data.fe_degree + 1,
2551 data.n_q_points_1d);
2552 else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2553 return Eval(data.shape_values,
2554 data.shape_gradients,
2555 data.shape_hessians,
2556 data.fe_degree + 1,
2557 data.n_q_points_1d);
2558 else
2559 {
2560 const unsigned int index =
2561 direction == 0 ? subface_index % 2 : subface_index / 2;
2562 return Eval(data.values_within_subface[index],
2563 data.gradients_within_subface[index],
2564 data.hessians_within_subface[index],
2565 data.fe_degree + 1,
2566 data.n_q_points_1d);
2567 }
2568 }
2569
2570 static void
2571 evaluate_in_face(
2572 const unsigned int n_components,
2573 const EvaluationFlags::EvaluationFlags evaluation_flag,
2574 const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
2575 Number * values_dofs,
2576 Number * values_quad,
2577 Number * gradients_quad,
2578 Number * hessians_quad,
2579 Number * scratch_data,
2580 const unsigned int subface_index)
2581 {
2582 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2583 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2584
2585 const std::size_t n_dofs = fe_degree > -1 ?
2586 Utilities::pow(fe_degree + 1, dim - 1) :
2587 Utilities::pow(data.fe_degree + 1, dim - 1);
2588 const std::size_t n_q_points =
2589 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2590 Utilities::pow(data.n_q_points_1d, dim - 1);
2591
2592 // keep a copy of the original pointer for the case of the Hessians
2593 Number *values_dofs_ptr = values_dofs;
2594
2595 if ((evaluation_flag & EvaluationFlags::values) != 0u &&
2596 ((evaluation_flag & EvaluationFlags::gradients) == 0u))
2597 for (unsigned int c = 0; c < n_components; ++c)
2598 {
2599 switch (dim)
2600 {
2601 case 3:
2602 eval0.template values<0, true, false>(values_dofs,
2603 values_quad);
2604 eval1.template values<1, true, false>(values_quad,
2605 values_quad);
2606 break;
2607 case 2:
2608 eval0.template values<0, true, false>(values_dofs,
2609 values_quad);
2610 break;
2611 case 1:
2612 values_quad[0] = values_dofs[0];
2613 break;
2614 default:
2615 Assert(false, ExcNotImplemented());
2616 }
2617 // Note: we always keep storage of values, 1st and 2nd derivatives
2618 // in an array
2619 values_dofs += 3 * n_dofs;
2620 values_quad += n_q_points;
2621 }
2622 else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2623 for (unsigned int c = 0; c < n_components; ++c)
2624 {
2625 switch (dim)
2626 {
2627 case 3:
2628 if (symmetric_evaluate &&
2629 use_collocation_evaluation(fe_degree, n_q_points_1d))
2630 {
2631 eval0.template values<0, true, false>(values_dofs,
2632 values_quad);
2633 eval0.template values<1, true, false>(values_quad,
2634 values_quad);
2635 EvaluatorTensorProduct<evaluate_evenodd,
2636 dim - 1,
2637 n_q_points_1d,
2638 n_q_points_1d,
2639 Number>
2640 eval_grad(AlignedVector<Number>(),
2641 data.shape_gradients_collocation_eo,
2642 AlignedVector<Number>());
2643 eval_grad.template gradients<0, true, false>(
2644 values_quad, gradients_quad);
2645 eval_grad.template gradients<1, true, false>(
2646 values_quad, gradients_quad + n_q_points);
2647 }
2648 else
2649 {
2650 // grad x
2651 eval0.template gradients<0, true, false>(values_dofs,
2652 scratch_data);
2653 eval1.template values<1, true, false>(scratch_data,
2654 gradients_quad);
2655
2656 // grad y
2657 eval0.template values<0, true, false>(values_dofs,
2658 scratch_data);
2659 eval1.template gradients<1, true, false>(scratch_data,
2660 gradients_quad +
2661 n_q_points);
2662
2663 if ((evaluation_flag & EvaluationFlags::values) != 0u)
2664 eval1.template values<1, true, false>(scratch_data,
2665 values_quad);
2666 }
2667 // grad z
2668 eval0.template values<0, true, false>(values_dofs + n_dofs,
2669 scratch_data);
2670 eval1.template values<1, true, false>(
2671 scratch_data, gradients_quad + (dim - 1) * n_q_points);
2672
2673 break;
2674 case 2:
2675 eval0.template values<0, true, false>(values_dofs + n_dofs,
2676 gradients_quad +
2677 n_q_points);
2678 eval0.template gradients<0, true, false>(values_dofs,
2679 gradients_quad);
2680 if ((evaluation_flag & EvaluationFlags::values) != 0u)
2681 eval0.template values<0, true, false>(values_dofs,
2682 values_quad);
2683 break;
2684 case 1:
2685 values_quad[0] = values_dofs[0];
2686 gradients_quad[0] = values_dofs[1];
2687 break;
2688 default:
2689 AssertThrow(false, ExcNotImplemented());
2690 }
2691 values_dofs += 3 * n_dofs;
2692 values_quad += n_q_points;
2693 gradients_quad += dim * n_q_points;
2694 }
2695
2696 if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
2697 {
2698 values_dofs = values_dofs_ptr;
2699 for (unsigned int c = 0; c < n_components; ++c)
2700 {
2701 switch (dim)
2702 {
2703 case 3:
2704 // grad xx
2705 eval0.template hessians<0, true, false>(values_dofs,
2706 scratch_data);
2707 eval1.template values<1, true, false>(scratch_data,
2708 hessians_quad);
2709
2710 // grad yy
2711 eval0.template values<0, true, false>(values_dofs,
2712 scratch_data);
2713 eval1.template hessians<1, true, false>(scratch_data,
2714 hessians_quad +
2715 n_q_points);
2716
2717 // grad zz
2718 eval0.template values<0, true, false>(values_dofs +
2719 2 * n_dofs,
2720 scratch_data);
2721 eval1.template values<1, true, false>(scratch_data,
2722 hessians_quad +
2723 2 * n_q_points);
2724
2725 // grad xy
2726 eval0.template gradients<0, true, false>(values_dofs,
2727 scratch_data);
2728 eval1.template gradients<1, true, false>(scratch_data,
2729 hessians_quad +
2730 3 * n_q_points);
2731
2732 // grad xz
2733 eval0.template gradients<0, true, false>(values_dofs +
2734 n_dofs,
2735 scratch_data);
2736 eval1.template values<1, true, false>(scratch_data,
2737 hessians_quad +
2738 4 * n_q_points);
2739
2740 // grad yz
2741 eval0.template values<0, true, false>(values_dofs + n_dofs,
2742 scratch_data);
2743 eval1.template gradients<1, true, false>(scratch_data,
2744 hessians_quad +
2745 5 * n_q_points);
2746
2747 break;
2748 case 2:
2749 // grad xx
2750 eval0.template hessians<0, true, false>(values_dofs,
2751 hessians_quad);
2752 // grad yy
2753 eval0.template values<0, true, false>(
2754 values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
2755 // grad xy
2756 eval0.template gradients<0, true, false>(
2757 values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
2758 break;
2759 case 1:
2760 hessians_quad[0] = values_dofs[2];
2761 break;
2762 default:
2763 AssertThrow(false, ExcNotImplemented());
2764 }
2765 values_dofs += 3 * n_dofs;
2766 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2767 }
2768 }
2769 }
2770
2771 static void
2772 integrate_in_face(
2773 const unsigned int n_components,
2774 const EvaluationFlags::EvaluationFlags integration_flag,
2775 const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
2776 Number * values_dofs,
2777 Number * values_quad,
2778 Number * gradients_quad,
2779 Number * hessians_quad,
2780 Number * scratch_data,
2781 const unsigned int subface_index)
2782 {
2783 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2784 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2785
2786 const std::size_t n_dofs =
2787 fe_degree > -1 ?
2788 Utilities::pow(fe_degree + 1, dim - 1) :
2789 (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
2790 const std::size_t n_q_points =
2791 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2792 Utilities::pow(data.n_q_points_1d, dim - 1);
2793
2794 // keep a copy of the original pointer for the case of the Hessians
2795 Number *values_dofs_ptr = values_dofs;
2796
2797 if ((integration_flag & EvaluationFlags::values) != 0u &&
2798 (integration_flag & EvaluationFlags::gradients) == 0u)
2799 for (unsigned int c = 0; c < n_components; ++c)
2800 {
2801 switch (dim)
2802 {
2803 case 3:
2804 eval1.template values<1, false, false>(values_quad,
2805 values_quad);
2806 eval0.template values<0, false, false>(values_quad,
2807 values_dofs);
2808 break;
2809 case 2:
2810 eval0.template values<0, false, false>(values_quad,
2811 values_dofs);
2812 break;
2813 case 1:
2814 values_dofs[0] = values_quad[0];
2815 break;
2816 default:
2817 Assert(false, ExcNotImplemented());
2818 }
2819 values_dofs += 3 * n_dofs;
2820 values_quad += n_q_points;
2821 }
2822 else if ((integration_flag & EvaluationFlags::gradients) != 0u)
2823 for (unsigned int c = 0; c < n_components; ++c)
2824 {
2825 switch (dim)
2826 {
2827 case 3:
2828 // grad z
2829 eval1.template values<1, false, false>(gradients_quad +
2830 2 * n_q_points,
2831 gradients_quad +
2832 2 * n_q_points);
2833 eval0.template values<0, false, false>(gradients_quad +
2834 2 * n_q_points,
2835 values_dofs + n_dofs);
2836 if (symmetric_evaluate &&
2837 use_collocation_evaluation(fe_degree, n_q_points_1d))
2838 {
2839 EvaluatorTensorProduct<evaluate_evenodd,
2840 dim - 1,
2841 n_q_points_1d,
2842 n_q_points_1d,
2843 Number>
2844 eval_grad(AlignedVector<Number>(),
2845 data.shape_gradients_collocation_eo,
2846 AlignedVector<Number>());
2847 if ((integration_flag & EvaluationFlags::values) != 0u)
2848 eval_grad.template gradients<1, false, true>(
2849 gradients_quad + n_q_points, values_quad);
2850 else
2851 eval_grad.template gradients<1, false, false>(
2852 gradients_quad + n_q_points, values_quad);
2853 eval_grad.template gradients<0, false, true>(
2854 gradients_quad, values_quad);
2855 eval0.template values<1, false, false>(values_quad,
2856 values_quad);
2857 eval0.template values<0, false, false>(values_quad,
2858 values_dofs);
2859 }
2860 else
2861 {
2862 if ((integration_flag & EvaluationFlags::values) != 0u)
2863 {
2864 eval1.template values<1, false, false>(values_quad,
2865 scratch_data);
2866 eval1.template gradients<1, false, true>(
2867 gradients_quad + n_q_points, scratch_data);
2868 }
2869 else
2870 eval1.template gradients<1, false, false>(
2871 gradients_quad + n_q_points, scratch_data);
2872
2873 // grad y
2874 eval0.template values<0, false, false>(scratch_data,
2875 values_dofs);
2876
2877 // grad x
2878 eval1.template values<1, false, false>(gradients_quad,
2879 scratch_data);
2880 eval0.template gradients<0, false, true>(scratch_data,
2881 values_dofs);
2882 }
2883 break;
2884 case 2:
2885 eval0.template values<0, false, false>(gradients_quad +
2886 n_q_points,
2887 values_dofs + n_dofs);
2888 eval0.template gradients<0, false, false>(gradients_quad,
2889 values_dofs);
2890 if ((integration_flag & EvaluationFlags::values) != 0u)
2891 eval0.template values<0, false, true>(values_quad,
2892 values_dofs);
2893 break;
2894 case 1:
2895 values_dofs[0] = values_quad[0];
2896 values_dofs[1] = gradients_quad[0];
2897 break;
2898 default:
2899 AssertThrow(false, ExcNotImplemented());
2900 }
2901 values_dofs += 3 * n_dofs;
2902 values_quad += n_q_points;
2903 gradients_quad += dim * n_q_points;
2904 }
2905
2906 if ((integration_flag & EvaluationFlags::hessians) != 0u)
2907 {
2908 values_dofs = values_dofs_ptr;
2909 for (unsigned int c = 0; c < n_components; ++c)
2910 {
2911 switch (dim)
2912 {
2913 case 3:
2914 // grad xx
2915 eval1.template values<1, false, false>(hessians_quad,
2916 scratch_data);
2917 if ((integration_flag & (EvaluationFlags::values |
2918 EvaluationFlags::gradients)) != 0u)
2919 eval0.template hessians<0, false, true>(scratch_data,
2920 values_dofs);
2921 else
2922 eval0.template hessians<0, false, false>(scratch_data,
2923 values_dofs);
2924
2925 // grad yy
2926 eval1.template hessians<1, false, false>(hessians_quad +
2927 n_q_points,
2928 scratch_data);
2929 eval0.template values<0, false, true>(scratch_data,
2930 values_dofs);
2931
2932 // grad zz
2933 eval1.template values<1, false, false>(hessians_quad +
2934 2 * n_q_points,
2935 scratch_data);
2936 eval0.template values<0, false, false>(scratch_data,
2937 values_dofs +
2938 2 * n_dofs);
2939
2940 // grad xy
2941 eval1.template gradients<1, false, false>(hessians_quad +
2942 3 * n_q_points,
2943 scratch_data);
2944 eval0.template gradients<0, false, true>(scratch_data,
2945 values_dofs);
2946
2947 // grad xz
2948 eval1.template values<1, false, false>(hessians_quad +
2949 4 * n_q_points,
2950 scratch_data);
2951 if ((integration_flag & EvaluationFlags::gradients) != 0u)
2952 eval0.template gradients<0, false, true>(scratch_data,
2953 values_dofs +
2954 n_dofs);
2955 else
2956 eval0.template gradients<0, false, false>(scratch_data,
2957 values_dofs +
2958 n_dofs);
2959
2960 // grad yz
2961 eval1.template gradients<1, false, false>(hessians_quad +
2962 5 * n_q_points,
2963 scratch_data);
2964 eval0.template values<0, false, true>(scratch_data,
2965 values_dofs + n_dofs);
2966
2967 break;
2968 case 2:
2969 // grad xx
2970 if ((integration_flag & (EvaluationFlags::values |
2971 EvaluationFlags::gradients)) != 0u)
2972 eval0.template hessians<0, false, true>(hessians_quad,
2973 values_dofs);
2974 else
2975 eval0.template hessians<0, false, false>(hessians_quad,
2976 values_dofs);
2977
2978 // grad yy
2979 eval0.template values<0, false, false>(
2980 hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
2981 // grad xy
2982 if ((integration_flag & EvaluationFlags::gradients) != 0u)
2983 eval0.template gradients<0, false, true>(
2984 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2985 else
2986 eval0.template gradients<0, false, false>(
2987 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2988 break;
2989 case 1:
2990 values_dofs[2] = hessians_quad[0];
2991 if ((integration_flag & EvaluationFlags::values) == 0u)
2992 values_dofs[0] = 0;
2993 if ((integration_flag & EvaluationFlags::gradients) == 0u)
2994 values_dofs[1] = 0;
2995 break;
2996 default:
2997 AssertThrow(false, ExcNotImplemented());
2998 }
2999 values_dofs += 3 * n_dofs;
3000 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
3001 }
3002 }
3003 }
3004 };
3005
3006 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
3007 struct FEFaceEvaluationImplRaviartThomas
3008 {
3009 using EvalGeneral = EvaluatorTensorProduct<evaluate_general,
3010 dim - 1,
3011 fe_degree,
3012 n_q_points_1d,
3013 Number>;
3014 template <typename EvalType>
3015 static EvalType
3016 create_evaluator_tensor_product(
3017 const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
3018 const unsigned int subface_index,
3019 const unsigned int direction)
3020 {
3021 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3022 return EvalType(data.shape_values,
3023 data.shape_gradients,
3024 data.shape_hessians);
3025 else
3026 {
3027 const unsigned int index =
3028 direction == 0 ? subface_index % 2 : subface_index / 2;
3029 return EvalType(data.values_within_subface[index],
3030 data.gradients_within_subface[index],
3031 data.hessians_within_subface[index]);
3032 }
3033 }
3034
3035 template <bool integrate>
3036 static void
3037 evaluate_or_integrate_in_face(
3038 const EvaluationFlags::EvaluationFlags evaluation_flag,
3039 Number * values_dofs,
3040 FEEvaluationData<dim, Number, true> & fe_eval,
3041 Number * scratch_data,
3042 const unsigned int subface_index,
3043 const unsigned int face_no)
3044 {
3045 const unsigned int face_direction = face_no / 2;
3046
3047 // We first evaluate the anisotropic faces, i.e the faces where
3048 // face_direction != component. Note that the call order here is not
3049 // important, since the pointers are shifted accordingly within the
3050 // function. However, this is the order in which the components will be in
3051 // the quadrature points. Furthermore, the isotropic faces have no "normal
3052 // direction" but we still pass in normal_dir = 2 since this is used for
3053 // the pointers.
3054 // -----------------------------------------------------------------------------------
3055 // | | Anisotropic faces | Isotropic faces|
3056 // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords |
3057 // | --------------------------------------------------------------------------------|
3058 // | 0 | 1, y, 0 | - | 0, y |
3059 // | 1 | 0, x, 0 | - | 1, x |
3060 // | --------------------------------------------------------------------------------|
3061 // | 0 | 1, yz, 0 | 2, yz, 1 | 0, yz |
3062 // | 1 | 2, zx, 0 | 0, zx, 1 | 1, zx |
3063 // | 2 | 0, xy, 0 | 1, xy, 1 | 2, xy |
3064 // -----------------------------------------------------------------------------------
3065 evaluate_in_face_apply<0>(values_dofs,
3066 fe_eval,
3067 scratch_data,
3068 evaluation_flag,
3069 face_direction,
3070 subface_index,
3071 std::integral_constant<bool, integrate>());
3072
3073 if (dim == 3)
3074 evaluate_in_face_apply<1>(values_dofs,
3075 fe_eval,
3076 scratch_data,
3077 evaluation_flag,
3078 face_direction,
3079 subface_index,
3080 std::integral_constant<bool, integrate>());
3081
3082 evaluate_in_face_apply<2>(values_dofs,
3083 fe_eval,
3084 scratch_data,
3085 evaluation_flag,
3086 face_direction,
3087 subface_index,
3088 std::integral_constant<bool, integrate>());
3089 }
3090
3091 /*
3092 * Helper function which applies the 1D kernels for on one
3093 * component in a face. normal_dir indicates the direction of the continuous
3094 * component of the RT space. std::integral_constant<bool, false> is the
3095 * evaluation path, and std::integral_constant<bool, true> below is the
3096 * integration path. These two functions can be fused together since all
3097 * offsets and pointers are the exact same.
3098 */
3099 template <int normal_dir>
3100 static inline void
3101 evaluate_in_face_apply(
3102 Number * values_dofs,
3103 FEEvaluationData<dim, Number, true> & fe_eval,
3104 Number * scratch_data,
3105 const EvaluationFlags::EvaluationFlags evaluation_flag,
3106 const unsigned int face_direction,
3107 const unsigned int subface_index,
3108 std::integral_constant<bool, false>)
3109 {
3110 using EvalNormal =
3111 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3112 dim - 1,
3113 (fe_degree == -1) ? 1 : fe_degree + 1,
3114 n_q_points_1d,
3115 Number,
3116 normal_dir>;
3117 using EvalTangent =
3118 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3119 dim - 1,
3120 (fe_degree == -1) ? 1 : fe_degree,
3121 n_q_points_1d,
3122 Number,
3123 normal_dir>;
3124
3125 using TempEval0 = typename std::
3126 conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3127 using TempEval1 = typename std::
3128 conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3129 using Eval0 = typename std::
3130 conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3131 using Eval1 = typename std::
3132 conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3133
3134 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3135 fe_eval.get_shape_info();
3136 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3137 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3138 subface_index,
3139 0);
3140 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3141 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3142 subface_index,
3143 1);
3144
3145 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3146 const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3147 const std::size_t n_dofs_normal =
3148 n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3149 const std::size_t dofs_stride =
3150 (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3151 n_dofs_tangent;
3152
3153 static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3154 {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3155 const unsigned int component =
3156 (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3157 0 :
3158 component_table[face_direction][normal_dir];
3159
3160 // Initial offsets
3161 values_dofs +=
3162 3 * ((component == 0) ?
3163 0 :
3164 ((component == 1) ?
3165 ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3166 ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3167 n_dofs_normal + n_dofs_tangent)));
3168 const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3169 Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3170 Number *gradients_quad =
3171 fe_eval.begin_gradients() + dim * n_q_points * shift;
3172 Number *hessians_quad =
3173 fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3174
3175 // Evaluation path
3176 if ((evaluation_flag & EvaluationFlags::values) &&
3177 !(evaluation_flag & EvaluationFlags::gradients))
3178 {
3179 switch (dim)
3180 {
3181 case 3:
3182 eval0.template values<0, true, false>(values_dofs, values_quad);
3183 eval1.template values<1, true, false>(values_quad, values_quad);
3184 break;
3185 case 2:
3186 eval0.template values<0, true, false>(values_dofs, values_quad);
3187 break;
3188 default:
3189 Assert(false, ExcNotImplemented());
3190 }
3191 }
3192 else if (evaluation_flag & EvaluationFlags::gradients)
3193 {
3194 switch (dim)
3195 {
3196 case 3:
3197 // grad x
3198 eval0.template gradients<0, true, false>(values_dofs,
3199 scratch_data);
3200 eval1.template values<1, true, false>(scratch_data,
3201 gradients_quad);
3202
3203 // grad y
3204 eval0.template values<0, true, false>(values_dofs,
3205 scratch_data);
3206 eval1.template gradients<1, true, false>(scratch_data,
3207 gradients_quad +
3208 n_q_points);
3209
3210 if (evaluation_flag & EvaluationFlags::values)
3211 eval1.template values<1, true, false>(scratch_data,
3212 values_quad);
3213
3214 // grad z
3215 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3216 scratch_data);
3217 eval1.template values<1, true, false>(scratch_data,
3218 gradients_quad +
3219 2 * n_q_points);
3220
3221 break;
3222 case 2:
3223 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3224 gradients_quad +
3225 n_q_points);
3226 eval0.template gradients<0, true, false>(values_dofs,
3227 gradients_quad);
3228 if ((evaluation_flag & EvaluationFlags::values))
3229 eval0.template values<0, true, false>(values_dofs,
3230 values_quad);
3231 break;
3232 default:
3233 AssertThrow(false, ExcNotImplemented());
3234 }
3235 }
3236
3237 if (evaluation_flag & EvaluationFlags::hessians)
3238 {
3239 switch (dim)
3240 {
3241 case 3:
3242 // grad xx
3243 eval0.template hessians<0, true, false>(values_dofs,
3244 scratch_data);
3245 eval1.template values<1, true, false>(scratch_data,
3246 hessians_quad);
3247
3248 // grad yy
3249 eval0.template values<0, true, false>(values_dofs,
3250 scratch_data);
3251 eval1.template hessians<1, true, false>(scratch_data,
3252 hessians_quad +
3253 n_q_points);
3254
3255 // grad zz
3256 eval0.template values<0, true, false>(values_dofs +
3257 2 * dofs_stride,
3258 scratch_data);
3259 eval1.template values<1, true, false>(scratch_data,
3260 hessians_quad +
3261 2 * n_q_points);
3262
3263 // grad xy
3264 eval0.template gradients<0, true, false>(values_dofs,
3265 scratch_data);
3266 eval1.template gradients<1, true, false>(scratch_data,
3267 hessians_quad +
3268 3 * n_q_points);
3269
3270 // grad xz
3271 eval0.template gradients<0, true, false>(values_dofs +
3272 dofs_stride,
3273 scratch_data);
3274 eval1.template values<1, true, false>(scratch_data,
3275 hessians_quad +
3276 4 * n_q_points);
3277
3278 // grad yz
3279 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3280 scratch_data);
3281 eval1.template gradients<1, true, false>(scratch_data,
3282 hessians_quad +
3283 5 * n_q_points);
3284
3285 break;
3286 case 2:
3287 // grad xx
3288 eval0.template hessians<0, true, false>(values_dofs,
3289 hessians_quad);
3290 // grad yy
3291 eval0.template values<0, true, false>(
3292 values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
3293 // grad xy
3294 eval0.template gradients<0, true, false>(
3295 values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
3296 break;
3297 default:
3298 AssertThrow(false, ExcNotImplemented());
3299 }
3300 }
3301 }
3302
3303 template <int normal_dir>
3304 static inline void
3305 evaluate_in_face_apply(
3306 Number * values_dofs,
3307 FEEvaluationData<dim, Number, true> & fe_eval,
3308 Number * scratch_data,
3309 const EvaluationFlags::EvaluationFlags evaluation_flag,
3310 const unsigned int face_direction,
3311 const unsigned int subface_index,
3312 std::integral_constant<bool, true>)
3313 {
3314 using EvalNormal =
3315 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3316 dim - 1,
3317 (fe_degree == -1) ? 1 : fe_degree + 1,
3318 n_q_points_1d,
3319 Number,
3320 normal_dir>;
3321 using EvalTangent =
3322 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3323 dim - 1,
3324 (fe_degree == -1) ? 1 : fe_degree,
3325 n_q_points_1d,
3326 Number,
3327 normal_dir>;
3328
3329 using TempEval0 = typename std::
3330 conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3331 using TempEval1 = typename std::
3332 conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3333 using Eval0 = typename std::
3334 conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3335 using Eval1 = typename std::
3336 conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3337
3338 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3339 fe_eval.get_shape_info();
3340 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3341 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3342 subface_index,
3343 0);
3344 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3345 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3346 subface_index,
3347 1);
3348
3349 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3350 const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3351 const std::size_t n_dofs_normal =
3352 n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3353 const std::size_t dofs_stride =
3354 (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3355 n_dofs_tangent;
3356
3357 static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3358 {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3359 const unsigned int component =
3360 (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3361 0 :
3362 component_table[face_direction][normal_dir];
3363
3364 // Initial offsets
3365 values_dofs +=
3366 3 * ((component == 0) ?
3367 0 :
3368 ((component == 1) ?
3369 ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3370 ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3371 n_dofs_normal + n_dofs_tangent)));
3372 const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3373 Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3374 Number *gradients_quad =
3375 fe_eval.begin_gradients() + dim * n_q_points * shift;
3376 Number *hessians_quad =
3377 fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3378
3379 // Integration path
3380 if ((evaluation_flag & EvaluationFlags::values) &&
3381 !(evaluation_flag & EvaluationFlags::gradients))
3382 {
3383 switch (dim)
3384 {
3385 case 3:
3386 eval1.template values<1, false, false>(values_quad,
3387 values_quad);
3388 eval0.template values<0, false, false>(values_quad,
3389 values_dofs);
3390 break;
3391 case 2:
3392 eval0.template values<0, false, false>(values_quad,
3393 values_dofs);
3394 break;
3395 default:
3396 Assert(false, ExcNotImplemented());
3397 }
3398 }
3399 else if (evaluation_flag & EvaluationFlags::gradients)
3400 {
3401 switch (dim)
3402 {
3403 case 3:
3404 // grad z
3405 eval1.template values<1, false, false>(gradients_quad +
3406 2 * n_q_points,
3407 gradients_quad +
3408 2 * n_q_points);
3409 eval0.template values<0, false, false>(
3410 gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
3411
3412 if (evaluation_flag & EvaluationFlags::values)
3413 {
3414 eval1.template values<1, false, false>(values_quad,
3415 scratch_data);
3416 eval1.template gradients<1, false, true>(gradients_quad +
3417 n_q_points,
3418 scratch_data);
3419 }
3420 else
3421 eval1.template gradients<1, false, false>(gradients_quad +
3422 n_q_points,
3423 scratch_data);
3424
3425 // grad y
3426 eval0.template values<0, false, false>(scratch_data,
3427 values_dofs);
3428
3429 // grad x
3430 eval1.template values<1, false, false>(gradients_quad,
3431 scratch_data);
3432 eval0.template gradients<0, false, true>(scratch_data,
3433 values_dofs);
3434
3435 break;
3436 case 2:
3437 eval0.template values<0, false, false>(
3438 gradients_quad + n_q_points, values_dofs + dofs_stride);
3439 eval0.template gradients<0, false, false>(gradients_quad,
3440 values_dofs);
3441 if (evaluation_flag & EvaluationFlags::values)
3442 eval0.template values<0, false, true>(values_quad,
3443 values_dofs);
3444 break;
3445 default:
3446 AssertThrow(false, ExcNotImplemented());
3447 }
3448 }
3449
3450 if (evaluation_flag & EvaluationFlags::hessians)
3451 {
3452 switch (dim)
3453 {
3454 case 3:
3455 // grad xx
3456 eval1.template values<1, false, false>(hessians_quad,
3457 scratch_data);
3458 if ((evaluation_flag &
3459 (EvaluationFlags::values | EvaluationFlags::gradients)))
3460 eval0.template hessians<0, false, true>(scratch_data,
3461 values_dofs);
3462 else
3463 eval0.template hessians<0, false, false>(scratch_data,
3464 values_dofs);
3465
3466 // grad yy
3467 eval1.template hessians<1, false, false>(hessians_quad +
3468 n_q_points,
3469 scratch_data);
3470 eval0.template values<0, false, true>(scratch_data,
3471 values_dofs);
3472
3473 // grad zz
3474 eval1.template values<1, false, false>(hessians_quad +
3475 2 * n_q_points,
3476 scratch_data);
3477 eval0.template values<0, false, false>(scratch_data,
3478 values_dofs +
3479 2 * dofs_stride);
3480
3481 // grad xy
3482 eval1.template gradients<1, false, false>(hessians_quad +
3483 3 * n_q_points,
3484 scratch_data);
3485 eval0.template gradients<0, false, true>(scratch_data,
3486 values_dofs);
3487
3488 // grad xz
3489 eval1.template values<1, false, false>(hessians_quad +
3490 4 * n_q_points,
3491 scratch_data);
3492 if ((evaluation_flag & EvaluationFlags::gradients))
3493 eval0.template gradients<0, false, true>(scratch_data,
3494 values_dofs +
3495 dofs_stride);
3496 else
3497 eval0.template gradients<0, false, false>(scratch_data,
3498 values_dofs +
3499 dofs_stride);
3500
3501 // grad yz
3502 eval1.template gradients<1, false, false>(hessians_quad +
3503 5 * n_q_points,
3504 scratch_data);
3505 eval0.template values<0, false, true>(scratch_data,
3506 values_dofs +
3507 dofs_stride);
3508
3509 break;
3510 case 2:
3511 // grad xx
3512 if (evaluation_flag &
3513 (EvaluationFlags::values | EvaluationFlags::gradients))
3514 eval0.template hessians<0, false, true>(hessians_quad,
3515 values_dofs);
3516 else
3517 eval0.template hessians<0, false, false>(hessians_quad,
3518 values_dofs);
3519
3520 // grad yy
3521 eval0.template values<0, false, false>(
3522 hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
3523 // grad xy
3524 if ((evaluation_flag & EvaluationFlags::gradients))
3525 eval0.template gradients<0, false, true>(
3526 hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3527 else
3528 eval0.template gradients<0, false, false>(
3529 hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3530 break;
3531 default:
3532 AssertThrow(false, ExcNotImplemented());
3533 }
3534 }
3535 }
3536 };
3537
3538
3539 template <int dim, int fe_degree, typename Number>
3540 struct FEFaceNormalEvaluationImpl
3541 {
3542 template <bool do_evaluate, bool add_into_output>
3543 static void
3544 interpolate(const unsigned int n_components,
3545 const EvaluationFlags::EvaluationFlags flags,
3546 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3547 const Number * input,
3548 Number * output,
3549 const unsigned int face_no)
3550 {
3551 Assert(static_cast<unsigned int>(fe_degree) ==
3552 shape_info.data.front().fe_degree ||
3553 fe_degree == -1,
3554 ExcInternalError());
3555 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
3556 interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
3557 n_components, input, output, flags, face_no, shape_info);
3558 else
3559 interpolate_generic<do_evaluate, add_into_output>(
3560 n_components,
3561 input,
3562 output,
3563 flags,
3564 face_no,
3565 shape_info.data.front().fe_degree + 1,
3566 shape_info.data.front().shape_data_on_face,
3567 shape_info.dofs_per_component_on_cell,
3568 3 * shape_info.dofs_per_component_on_face);
3569 }
3570
3574 template <bool do_evaluate, bool add_into_output>
3575 static void
3576 interpolate_quadrature(
3577 const unsigned int n_components,
3578 const EvaluationFlags::EvaluationFlags flags,
3579 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3580 const Number * input,
3581 Number * output,
3582 const unsigned int face_no)
3583 {
3584 Assert(static_cast<unsigned int>(fe_degree + 1) ==
3585 shape_info.data.front().n_q_points_1d ||
3586 fe_degree == -1,
3587 ExcInternalError());
3588
3589 interpolate_generic<do_evaluate, add_into_output>(
3590 n_components,
3591 input,
3592 output,
3593 flags,
3594 face_no,
3595 shape_info.data.front().quadrature.size(),
3596 shape_info.data.front().quadrature_data_on_face,
3597 shape_info.n_q_points,
3598 shape_info.n_q_points_face);
3599 }
3600
3601 private:
3602 template <bool do_evaluate, bool add_into_output, int face_direction = 0>
3603 static void
3604 interpolate_generic(const unsigned int n_components,
3605 const Number * input,
3606 Number * output,
3607 const EvaluationFlags::EvaluationFlags flag,
3608 const unsigned int face_no,
3609 const unsigned int n_points_1d,
3610 const std::array<AlignedVector<Number>, 2> &shape_data,
3611 const unsigned int dofs_per_component_on_cell,
3612 const unsigned int dofs_per_component_on_face)
3613 {
3614 if (face_direction == face_no / 2)
3615 {
3616 EvaluatorTensorProduct<evaluate_general,
3617 dim,
3618 fe_degree + 1,
3619 0,
3620 Number>
3621 evalf(shape_data[face_no % 2],
3622 AlignedVector<Number>(),
3623 AlignedVector<Number>(),
3624 n_points_1d,
3625 0);
3626
3627 const unsigned int in_stride = do_evaluate ?
3628 dofs_per_component_on_cell :
3629 dofs_per_component_on_face;
3630 const unsigned int out_stride = do_evaluate ?
3631 dofs_per_component_on_face :
3632 dofs_per_component_on_cell;
3633
3634 for (unsigned int c = 0; c < n_components; ++c)
3635 {
3636 if (flag & EvaluationFlags::hessians)
3637 evalf.template apply_face<face_direction,
3638 do_evaluate,
3639 add_into_output,
3640 2>(input, output);
3641 else if (flag & EvaluationFlags::gradients)
3642 evalf.template apply_face<face_direction,
3643 do_evaluate,
3644 add_into_output,
3645 1>(input, output);
3646 else
3647 evalf.template apply_face<face_direction,
3648 do_evaluate,
3649 add_into_output,
3650 0>(input, output);
3651 input += in_stride;
3652 output += out_stride;
3653 }
3654 }
3655 else if (face_direction < dim)
3656 {
3657 interpolate_generic<do_evaluate,
3658 add_into_output,
3659 std::min(face_direction + 1, dim - 1)>(
3660 n_components,
3661 input,
3662 output,
3663 flag,
3664 face_no,
3665 n_points_1d,
3666 shape_data,
3667 dofs_per_component_on_cell,
3668 dofs_per_component_on_face);
3669 }
3670 }
3671
3672 template <typename EvalType>
3673 static EvalType
3674 create_evaluator_tensor_product(
3675 const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
3676 const unsigned int face_no)
3677 {
3678 return EvalType(data.shape_data_on_face[face_no % 2],
3679 AlignedVector<Number>(),
3680 AlignedVector<Number>());
3681 }
3682
3683 template <bool do_evaluate,
3684 bool add_into_output,
3685 int face_direction = 0,
3686 int max_derivative = 0>
3687 static void
3688 interpolate_generic_raviart_thomas(
3689 const unsigned int n_components,
3690 const Number * input,
3691 Number * output,
3692 const EvaluationFlags::EvaluationFlags flag,
3693 const unsigned int face_no,
3694 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
3695 {
3696 if (dim == 1)
3697 {
3698 // This should never happen since the FE_RaviartThomasNodal is not
3699 // defined for dim = 1. It prevents compiler warnings of infinite
3700 // recursion.
3701 Assert(false, ExcInternalError());
3702 return;
3703 }
3704
3705 bool increase_max_der = false;
3706 if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
3707 (flag & EvaluationFlags::gradients && max_derivative < 1))
3708 increase_max_der = true;
3709
3710 if (face_direction == face_no / 2 && !increase_max_der)
3711 {
3712 interpolate_generic_raviart_thomas_apply_face<do_evaluate,
3713 add_into_output,
3714 face_direction,
3715 max_derivative>(
3716 shape_info, face_no, input, output);
3717 }
3718 else if (face_direction == face_no / 2)
3719 {
3720 // Only increase max_derivative
3721 interpolate_generic_raviart_thomas<do_evaluate,
3722 add_into_output,
3723 face_direction,
3724 std::min(max_derivative + 1, 2)>(
3725 n_components, input, output, flag, face_no, shape_info);
3726 }
3727 else if (face_direction < dim)
3728 {
3729 if (increase_max_der)
3730 {
3731 interpolate_generic_raviart_thomas<
3732 do_evaluate,
3733 add_into_output,
3734 std::min(face_direction + 1, dim - 1),
3735 std::min(max_derivative + 1, 2)>(
3736 n_components, input, output, flag, face_no, shape_info);
3737 }
3738 else
3739 {
3740 interpolate_generic_raviart_thomas<do_evaluate,
3741 add_into_output,
3742 std::min(face_direction + 1,
3743 dim - 1),
3744 max_derivative>(
3745 n_components, input, output, flag, face_no, shape_info);
3746 }
3747 }
3748 }
3749
3750 /* Help function for interpolate_generic_raviart_thomas */
3751 template <bool do_evaluate,
3752 bool add_into_output,
3753 int face_direction,
3754 int max_derivative>
3755 static inline void
3756 interpolate_generic_raviart_thomas_apply_face(
3757 const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3758 const unsigned int face_no,
3759 const Number * input,
3760 Number * output)
3761 {
3762 // These types are evaluators in either normal or tangential direction
3763 // depending on the face direction, with different normal directions for
3764 // the different components.
3765 using Evalf0 = typename std::conditional<
3766 face_direction == 0,
3767 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3768 dim,
3769 (fe_degree == -1) ? 1 : fe_degree + 1,
3770 0,
3771 Number,
3772 0>,
3773 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3774 dim,
3775 (fe_degree == -1) ? 1 : fe_degree,
3776 0,
3777 Number,
3778 0>>::type;
3779 using Evalf1 = typename std::conditional<
3780 face_direction == 1,
3781 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3782 dim,
3783 (fe_degree == -1) ? 1 : fe_degree + 1,
3784 0,
3785 Number,
3786 1>,
3787 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3788 dim,
3789 (fe_degree == -1) ? 1 : fe_degree,
3790 0,
3791 Number,
3792 1>>::type;
3793 using Evalf2 = typename std::conditional<
3794 face_direction == 2,
3795 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3796 dim,
3797 (fe_degree == -1) ? 1 : fe_degree + 1,
3798 0,
3799 Number,
3800 2>,
3801 EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
3802 dim,
3803 (fe_degree == -1) ? 1 : fe_degree,
3804 0,
3805 Number,
3806 2>>::type;
3807
3808 Evalf0 evalf0 =
3809 create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
3810 shape_info.data[0] :
3811 shape_info.data[1],
3812 face_no);
3813 Evalf1 evalf1 =
3814 create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
3815 shape_info.data[0] :
3816 shape_info.data[1],
3817 face_no);
3818 Evalf2 evalf2 =
3819 create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
3820 shape_info.data[0] :
3821 shape_info.data[1],
3822 face_no);
3823
3824 const unsigned int dofs_per_component_on_cell =
3825 shape_info.dofs_per_component_on_cell;
3826 const unsigned int dofs_per_component_on_face =
3827 3 * shape_info.dofs_per_component_on_face;
3828
3829 // NOTE! dofs_per_component_on_face is in the tangent direction,
3830 // i.e (fe.degree+1)*fe.degree. Normal faces are only
3831 // fe.degree*fe.degree
3832 const unsigned int in_stride =
3833 do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
3834 const unsigned int out_stride =
3835 do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
3836
3837 const unsigned int in_stride_after_normal =
3838 do_evaluate ?
3839 dofs_per_component_on_cell :
3840 dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
3841 const unsigned int out_stride_after_normal =
3842 do_evaluate ?
3843 dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
3844 dofs_per_component_on_cell;
3845
3846 evalf0.template apply_face<face_direction,
3847 do_evaluate,
3848 add_into_output,
3849 max_derivative>(input, output);
3850 // stride to next component
3851 input += (face_direction == 0) ? in_stride_after_normal : in_stride;
3852 output += (face_direction == 0) ? out_stride_after_normal : out_stride;
3853
3854 evalf1.template apply_face<face_direction,
3855 do_evaluate,
3856 add_into_output,
3857 max_derivative>(input, output);
3858
3859 if (dim == 3)
3860 {
3861 // stride to next component
3862 input += (face_direction == 1) ? in_stride_after_normal : in_stride;
3863 output +=
3864 (face_direction == 1) ? out_stride_after_normal : out_stride;
3865
3866 evalf2.template apply_face<face_direction,
3867 do_evaluate,
3868 add_into_output,
3869 max_derivative>(input, output);
3870 }
3871 }
3872 };
3873
3874
3875
3876 // internal helper function for reading data; base version of different types
3877 template <typename VectorizedArrayType, typename Number2>
3878 void
3879 do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
3880 {
3881 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3882 dst[v] = src_ptr[v];
3883 }
3884
3885
3886
3887 // internal helper function for reading data; specialized version where we
3888 // can use a dedicated load function
3889 template <typename Number, std::size_t width>
3890 void
3891 do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
3892 {
3893 dst.load(src_ptr);
3894 }
3895
3896
3897
3898 // internal helper function for reading data; base version of different types
3899 template <typename VectorizedArrayType, typename Number2>
3900 void
3901 do_vectorized_gather(const Number2 * src_ptr,
3902 const unsigned int * indices,
3903 VectorizedArrayType &dst)
3904 {
3905 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3906 dst[v] = src_ptr[indices[v]];
3907 }
3908
3909
3910
3911 // internal helper function for reading data; specialized version where we
3912 // can use a dedicated gather function
3913 template <typename Number, std::size_t width>
3914 void
3915 do_vectorized_gather(const Number * src_ptr,
3916 const unsigned int * indices,
3917 VectorizedArray<Number, width> &dst)
3918 {
3919 dst.gather(src_ptr, indices);
3920 }
3921
3922
3923
3924 // internal helper function for reading data; base version of different types
3925 template <typename VectorizedArrayType, typename Number2>
3926 void
3927 do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
3928 {
3929 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3930 dst_ptr[v] += src[v];
3931 }
3932
3933
3934
3935 // internal helper function for reading data; specialized version where we
3936 // can use a dedicated load function
3937 template <typename Number, std::size_t width>
3938 void
3939 do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
3940 {
3941 VectorizedArray<Number, width> tmp;
3942 tmp.load(dst_ptr);
3943 (tmp + src).store(dst_ptr);
3944 }
3945
3946
3947
3948 // internal helper function for reading data; base version of different types
3949 template <typename VectorizedArrayType, typename Number2>
3950 void
3951 do_vectorized_scatter_add(const VectorizedArrayType src,
3952 const unsigned int * indices,
3953 Number2 * dst_ptr)
3954 {
3955 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3956 dst_ptr[indices[v]] += src[v];
3957 }
3958
3959
3960
3961 // internal helper function for reading data; specialized version where we
3962 // can use a dedicated gather function
3963 template <typename Number, std::size_t width>
3964 void
3965 do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
3966 const unsigned int * indices,
3967 Number * dst_ptr)
3968 {
3969#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
3970 for (unsigned int v = 0; v < width; ++v)
3971 dst_ptr[indices[v]] += src[v];
3972#else
3973 VectorizedArray<Number, width> tmp;
3974 tmp.gather(dst_ptr, indices);
3975 (tmp + src).scatter(indices, dst_ptr);
3976#endif
3977 }
3978
3979
3980
3981 template <typename Number>
3982 void
3983 adjust_for_face_orientation(const unsigned int dim,
3984 const unsigned int n_components,
3985 const EvaluationFlags::EvaluationFlags flag,
3986 const unsigned int *orientation,
3987 const bool integrate,
3988 const std::size_t n_q_points,
3989 Number * tmp_values,
3990 Number * values_quad,
3991 Number * gradients_quad,
3992 Number * hessians_quad)
3993 {
3994 for (unsigned int c = 0; c < n_components; ++c)
3995 {
3996 if (flag & EvaluationFlags::values)
3997 {
3998 if (integrate)
3999 for (unsigned int q = 0; q < n_q_points; ++q)
4000 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
4001 else
4002 for (unsigned int q = 0; q < n_q_points; ++q)
4003 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
4004 for (unsigned int q = 0; q < n_q_points; ++q)
4005 values_quad[c * n_q_points + q] = tmp_values[q];
4006 }
4007 if (flag & EvaluationFlags::gradients)
4008 for (unsigned int d = 0; d < dim; ++d)
4009 {
4010 if (integrate)
4011 for (unsigned int q = 0; q < n_q_points; ++q)
4012 tmp_values[q] =
4013 gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
4014 else
4015 for (unsigned int q = 0; q < n_q_points; ++q)
4016 tmp_values[orientation[q]] =
4017 gradients_quad[(c * dim + d) * n_q_points + q];
4018 for (unsigned int q = 0; q < n_q_points; ++q)
4019 gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
4020 }
4021 if (flag & EvaluationFlags::hessians)
4022 {
4023 const unsigned int hdim = (dim * (dim + 1)) / 2;
4024 for (unsigned int d = 0; d < hdim; ++d)
4025 {
4026 if (integrate)
4027 for (unsigned int q = 0; q < n_q_points; ++q)
4028 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4029 orientation[q]];
4030 else
4031 for (unsigned int q = 0; q < n_q_points; ++q)
4032 tmp_values[orientation[q]] =
4033 hessians_quad[(c * hdim + d) * n_q_points + q];
4034 for (unsigned int q = 0; q < n_q_points; ++q)
4035 hessians_quad[(c * hdim + d) * n_q_points + q] =
4036 tmp_values[q];
4037 }
4038 }
4039 }
4040 }
4041
4042
4043
4044 template <typename Number, typename VectorizedArrayType>
4045 void
4046 adjust_for_face_orientation_per_lane(
4047 const unsigned int dim,
4048 const unsigned int n_components,
4049 const unsigned int v,
4050 const EvaluationFlags::EvaluationFlags flag,
4051 const unsigned int * orientation,
4052 const bool integrate,
4053 const std::size_t n_q_points,
4054 Number * tmp_values,
4055 VectorizedArrayType * values_quad,
4056 VectorizedArrayType * gradients_quad = nullptr,
4057 VectorizedArrayType * hessians_quad = nullptr)
4058 {
4059 for (unsigned int c = 0; c < n_components; ++c)
4060 {
4061 if (flag & EvaluationFlags::values)
4062 {
4063 if (integrate)
4064 for (unsigned int q = 0; q < n_q_points; ++q)
4065 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
4066 else
4067 for (unsigned int q = 0; q < n_q_points; ++q)
4068 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
4069 for (unsigned int q = 0; q < n_q_points; ++q)
4070 values_quad[c * n_q_points + q][v] = tmp_values[q];
4071 }
4072 if (flag & EvaluationFlags::gradients)
4073 for (unsigned int d = 0; d < dim; ++d)
4074 {
4075 Assert(gradients_quad != nullptr, ExcInternalError());
4076 if (integrate)
4077 for (unsigned int q = 0; q < n_q_points; ++q)
4078 tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
4079 orientation[q]][v];
4080 else
4081 for (unsigned int q = 0; q < n_q_points; ++q)
4082 tmp_values[orientation[q]] =
4083 gradients_quad[(c * dim + d) * n_q_points + q][v];
4084 for (unsigned int q = 0; q < n_q_points; ++q)
4085 gradients_quad[(c * dim + d) * n_q_points + q][v] =
4086 tmp_values[q];
4087 }
4088 if (flag & EvaluationFlags::hessians)
4089 {
4090 Assert(hessians_quad != nullptr, ExcInternalError());
4091 const unsigned int hdim = (dim * (dim + 1)) / 2;
4092 for (unsigned int d = 0; d < hdim; ++d)
4093 {
4094 if (integrate)
4095 for (unsigned int q = 0; q < n_q_points; ++q)
4096 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4097 orientation[q]][v];
4098 else
4099 for (unsigned int q = 0; q < n_q_points; ++q)
4100 tmp_values[orientation[q]] =
4101 hessians_quad[(c * hdim + d) * n_q_points + q][v];
4102 for (unsigned int q = 0; q < n_q_points; ++q)
4103 hessians_quad[(c * hdim + d) * n_q_points + q][v] =
4104 tmp_values[q];
4105 }
4106 }
4107 }
4108 }
4109
4110
4111
4112 template <int dim, typename Number>
4113 struct FEFaceEvaluationImplEvaluateSelector
4114 {
4115 template <int fe_degree, int n_q_points_1d>
4116 static bool
4117 run(const unsigned int n_components,
4118 const EvaluationFlags::EvaluationFlags evaluation_flag,
4119 const Number * values_dofs,
4120 FEEvaluationData<dim, Number, true> & fe_eval)
4121 {
4122 const auto &shape_info = fe_eval.get_shape_info();
4123 const auto &shape_data = shape_info.data.front();
4124
4125 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4126 {
4127 Assert((fe_eval.get_dof_access_index() ==
4128 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
4129 fe_eval.is_interior_face() == false) == false,
4130 ExcNotImplemented());
4131
4132 const unsigned int face_no = fe_eval.get_face_no();
4133 const unsigned int face_orientation = fe_eval.get_face_orientation();
4134 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4135 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4136
4137 using Eval =
4138 EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
4139
4140 if (evaluation_flag & EvaluationFlags::values)
4141 {
4142 const auto shape_values =
4143 &shape_data.shape_values_face(face_no, face_orientation, 0);
4144
4145 auto values_quad_ptr = fe_eval.begin_values();
4146 auto values_dofs_actual_ptr = values_dofs;
4147
4148 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4149 for (unsigned int c = 0; c < n_components; ++c)
4150 {
4151 eval.template values<0, true, false>(values_dofs_actual_ptr,
4152 values_quad_ptr);
4153
4154 values_quad_ptr += n_q_points;
4155 values_dofs_actual_ptr += n_dofs;
4156 }
4157 }
4158
4159 if (evaluation_flag & EvaluationFlags::gradients)
4160 {
4161 auto gradients_quad_ptr = fe_eval.begin_gradients();
4162 auto values_dofs_actual_ptr = values_dofs;
4163
4164 std::array<const Number *, dim> shape_gradients;
4165 for (unsigned int d = 0; d < dim; ++d)
4166 shape_gradients[d] = &shape_data.shape_gradients_face(
4167 face_no, face_orientation, d, 0);
4168
4169 for (unsigned int c = 0; c < n_components; ++c)
4170 {
4171 for (unsigned int d = 0; d < dim; ++d)
4172 {
4173 Eval eval(nullptr,
4174 shape_gradients[d],
4175 nullptr,
4176 n_dofs,
4177 n_q_points);
4178
4179 eval.template gradients<0, true, false>(
4180 values_dofs_actual_ptr, gradients_quad_ptr);
4181
4182 gradients_quad_ptr += n_q_points;
4183 }
4184 values_dofs_actual_ptr += n_dofs;
4185 }
4186 }
4187
4188 Assert(!(evaluation_flag & EvaluationFlags::hessians),
4189 ExcNotImplemented());
4190
4191 return true;
4192 }
4193
4194 const unsigned int dofs_per_face =
4195 fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4196 Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4197
4198 // Note: we always keep storage of values, 1st and 2nd derivatives in an
4199 // array, so reserve space for all three here
4200 Number *temp = fe_eval.get_scratch_data().begin();
4201 Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4202
4203 bool use_vectorization = true;
4204
4205 if (fe_eval.get_dof_access_index() ==
4206 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
4207 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4208 use_vectorization =
4209 fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
4210 std::all_of(fe_eval.get_cell_ids().begin() + 1,
4211 fe_eval.get_cell_ids().end(),
4212 [&](const auto &v) {
4213 return v == fe_eval.get_cell_ids()[0] ||
4214 v == numbers::invalid_unsigned_int;
4215 });
4216
4217 if (use_vectorization == false)
4218 {
4219 for (unsigned int v = 0; v < Number::size(); ++v)
4220 {
4221 // the loop breaks once an invalid_unsigned_int is hit for
4222 // all cases except the exterior faces in the ECL loop (where
4223 // some faces might be at the boundaries but others not)
4224 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4225 {
4226 for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4227 ++i)
4228 temp[i][v] = 0;
4229 continue;
4230 }
4231
4232 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
4233 template interpolate<true, false>(n_components,
4234 evaluation_flag,
4235 shape_info,
4236 values_dofs,
4237 scratch_data,
4238 fe_eval.get_face_no(v));
4239
4240 for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4241 ++i)
4242 temp[i][v] = scratch_data[i][v];
4243 }
4244 }
4245 else
4246 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
4247 template interpolate<true, false>(n_components,
4248 evaluation_flag,
4249 shape_info,
4250 values_dofs,
4251 temp,
4252 fe_eval.get_face_no());
4253
4254 const unsigned int subface_index = fe_eval.get_subface_index();
4255 constexpr unsigned int n_q_points_1d_actual =
4256 fe_degree > -1 ? n_q_points_1d : 0;
4257
4258 if (fe_degree >= 1 &&
4259 shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4260 {
4261 FEFaceEvaluationImplRaviartThomas<dim,
4262 (fe_degree == -1) ? 1 : fe_degree,
4263 (n_q_points_1d < 1) ? 1 :
4264 n_q_points_1d,
4265 Number>::
4266 template evaluate_or_integrate_in_face<false>(
4267 evaluation_flag,
4268 temp,
4269 fe_eval,
4270 scratch_data,
4271 subface_index,
4272 fe_eval.get_face_no());
4273 }
4274 else if (fe_degree > -1 &&
4275 subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4276 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4277 FEFaceEvaluationImpl<true,
4278 dim,
4279 fe_degree,
4280 n_q_points_1d_actual,
4281 Number>::evaluate_in_face(n_components,
4282 evaluation_flag,
4283 shape_data,
4284 temp,
4285 fe_eval.begin_values(),
4286 fe_eval
4287 .begin_gradients(),
4288 fe_eval.begin_hessians(),
4289 scratch_data,
4290 subface_index);
4291 else
4292 FEFaceEvaluationImpl<false,
4293 dim,
4294 fe_degree,
4295 n_q_points_1d_actual,
4296 Number>::evaluate_in_face(n_components,
4297 evaluation_flag,
4298 shape_data,
4299 temp,
4300 fe_eval.begin_values(),
4301 fe_eval
4302 .begin_gradients(),
4303 fe_eval.begin_hessians(),
4304 scratch_data,
4305 subface_index);
4306
4307 if (use_vectorization == false)
4308 {
4309 for (unsigned int v = 0; v < Number::size(); ++v)
4310 {
4311 // the loop breaks once an invalid_unsigned_int is hit for
4312 // all cases except the exterior faces in the ECL loop (where
4313 // some faces might be at the boundaries but others not)
4314 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4315 continue;
4316
4317 if (fe_eval.get_face_orientation(v) != 0)
4318 adjust_for_face_orientation_per_lane(
4319 dim,
4320 n_components,
4321 v,
4322 evaluation_flag,
4323 &fe_eval.get_shape_info().face_orientations_quad(
4324 fe_eval.get_face_orientation(v), 0),
4325 false,
4326 shape_info.n_q_points_face,
4327 &temp[0][0],
4328 fe_eval.begin_values(),
4329 fe_eval.begin_gradients(),
4330 fe_eval.begin_hessians());
4331 }
4332 }
4333 else if (fe_eval.get_face_orientation() != 0)
4334 adjust_for_face_orientation(
4335 dim,
4336 n_components,
4337 evaluation_flag,
4338 &fe_eval.get_shape_info().face_orientations_quad(
4339 fe_eval.get_face_orientation(), 0),
4340 false,
4341 shape_info.n_q_points_face,
4342 temp,
4343 fe_eval.begin_values(),
4344 fe_eval.begin_gradients(),
4345 fe_eval.begin_hessians());
4346
4347 return false;
4348 }
4349 };
4350
4351
4352
4353 template <int dim, typename Number>
4354 struct FEFaceEvaluationImplIntegrateSelector
4355 {
4356 template <int fe_degree, int n_q_points_1d>
4357 static bool
4358 run(const unsigned int n_components,
4359 const EvaluationFlags::EvaluationFlags integration_flag,
4360 Number * values_dofs,
4361 FEEvaluationData<dim, Number, true> & fe_eval)
4362 {
4363 const auto &shape_info = fe_eval.get_shape_info();
4364 const auto &shape_data = shape_info.data.front();
4365
4366 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4367 {
4368 Assert((fe_eval.get_dof_access_index() ==
4369 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
4370 fe_eval.is_interior_face() == false) == false,
4371 ExcNotImplemented());
4372
4373 const unsigned int face_no = fe_eval.get_face_no();
4374 const unsigned int face_orientation = fe_eval.get_face_orientation();
4375 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4376 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4377
4378 using Eval =
4379 EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
4380
4381 if (integration_flag & EvaluationFlags::values)
4382 {
4383 const auto shape_values =
4384 &shape_data.shape_values_face(face_no, face_orientation, 0);
4385
4386 auto values_quad_ptr = fe_eval.begin_values();
4387 auto values_dofs_actual_ptr = values_dofs;
4388
4389 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4390 for (unsigned int c = 0; c < n_components; ++c)
4391 {
4392 eval.template values<0, false, false>(values_quad_ptr,
4393 values_dofs_actual_ptr);
4394
4395 values_quad_ptr += n_q_points;
4396 values_dofs_actual_ptr += n_dofs;
4397 }
4398 }
4399
4400 if (integration_flag & EvaluationFlags::gradients)
4401 {
4402 auto gradients_quad_ptr = fe_eval.begin_gradients();
4403 auto values_dofs_actual_ptr = values_dofs;
4404
4405 std::array<const Number *, dim> shape_gradients;
4406 for (unsigned int d = 0; d < dim; ++d)
4407 shape_gradients[d] = &shape_data.shape_gradients_face(
4408 face_no, face_orientation, d, 0);
4409
4410 for (unsigned int c = 0; c < n_components; ++c)
4411 {
4412 for (unsigned int d = 0; d < dim; ++d)
4413 {
4414 Eval eval(nullptr,
4415 shape_gradients[d],
4416 nullptr,
4417 n_dofs,
4418 n_q_points);
4419
4420 if (!(integration_flag & EvaluationFlags::values) &&
4421 d == 0)
4422 eval.template gradients<0, false, false>(
4423 gradients_quad_ptr, values_dofs_actual_ptr);
4424 else
4425 eval.template gradients<0, false, true>(
4426 gradients_quad_ptr, values_dofs_actual_ptr);
4427
4428 gradients_quad_ptr += n_q_points;
4429 }
4430 values_dofs_actual_ptr += n_dofs;
4431 }
4432 }
4433
4434 Assert(!(integration_flag & EvaluationFlags::hessians),
4435 ExcNotImplemented());
4436
4437 return true;
4438 }
4439
4440 const unsigned int dofs_per_face =
4441 fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4442 Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4443
4444 Number *temp = fe_eval.get_scratch_data().begin();
4445 Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4446
4447 bool use_vectorization = true;
4448
4449 if (fe_eval.get_dof_access_index() ==
4450 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
4451 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4452 use_vectorization =
4453 fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
4454 std::all_of(fe_eval.get_cell_ids().begin() + 1,
4455 fe_eval.get_cell_ids().end(),
4456 [&](const auto &v) {
4457 return v == fe_eval.get_cell_ids()[0] ||
4458 v == numbers::invalid_unsigned_int;
4459 });
4460
4461 if (use_vectorization == false)
4462 {
4463 for (unsigned int v = 0; v < Number::size(); ++v)
4464 {
4465 // the loop breaks once an invalid_unsigned_int is hit for
4466 // all cases except the exterior faces in the ECL loop (where
4467 // some faces might be at the boundaries but others not)
4468 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4469 continue;
4470
4471 if (fe_eval.get_face_orientation(v) != 0)
4472 adjust_for_face_orientation_per_lane(
4473 dim,
4474 n_components,
4475 v,
4476 integration_flag,
4477 &fe_eval.get_shape_info().face_orientations_quad(
4478 fe_eval.get_face_orientation(v), 0),
4479 true,
4480 shape_info.n_q_points_face,
4481 &temp[0][0],
4482 fe_eval.begin_values(),
4483 fe_eval.begin_gradients(),
4484 fe_eval.begin_hessians());
4485 }
4486 }
4487 else if (fe_eval.get_face_orientation() != 0)
4488 adjust_for_face_orientation(
4489 dim,
4490 n_components,
4491 integration_flag,
4492 &fe_eval.get_shape_info().face_orientations_quad(
4493 fe_eval.get_face_orientation(), 0),
4494 true,
4495 shape_info.n_q_points_face,
4496 temp,
4497 fe_eval.begin_values(),
4498 fe_eval.begin_gradients(),
4499 fe_eval.begin_hessians());
4500
4501 const unsigned int n_q_points_1d_actual =
4502 fe_degree > -1 ? n_q_points_1d : 0;
4503 const unsigned int subface_index = fe_eval.get_subface_index();
4504
4505 if (fe_degree >= 1 &&
4506 shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4507 {
4508 FEFaceEvaluationImplRaviartThomas<dim,
4509 (fe_degree == -1) ? 1 : fe_degree,
4510 (n_q_points_1d < 1) ? 1 :
4511 n_q_points_1d,
4512 Number>::
4513 template evaluate_or_integrate_in_face<true>(integration_flag,
4514 temp,
4515 fe_eval,
4516 scratch_data,
4517 subface_index,
4518 fe_eval.get_face_no());
4519 }
4520 else if (fe_degree > -1 &&
4521 fe_eval.get_subface_index() >=
4522 GeometryInfo<dim - 1>::max_children_per_cell &&
4523 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4524 FEFaceEvaluationImpl<
4525 true,
4526 dim,
4527 fe_degree,
4528 n_q_points_1d_actual,
4529 Number>::integrate_in_face(n_components,
4530 integration_flag,
4531 shape_data,
4532 temp,
4533 fe_eval.begin_values(),
4534 fe_eval.begin_gradients(),
4535 fe_eval.begin_hessians(),
4536 scratch_data,
4537 subface_index);
4538 else
4539 FEFaceEvaluationImpl<
4540 false,
4541 dim,
4542 fe_degree,
4543 n_q_points_1d_actual,
4544 Number>::integrate_in_face(n_components,
4545 integration_flag,
4546 shape_data,
4547 temp,
4548 fe_eval.begin_values(),
4549 fe_eval.begin_gradients(),
4550 fe_eval.begin_hessians(),
4551 scratch_data,
4552 subface_index);
4553
4554 if (use_vectorization == false)
4555 {
4556 for (unsigned int v = 0; v < Number::size(); ++v)
4557 {
4558 // the loop breaks once an invalid_unsigned_int is hit for
4559 // all cases except the exterior faces in the ECL loop (where
4560 // some faces might be at the boundaries but others not)
4561 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4562 continue;
4563
4564 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
4565 template interpolate<false, false>(n_components,
4566 integration_flag,
4567 shape_info,
4568 values_dofs,
4569 scratch_data,
4570 fe_eval.get_face_no(v));
4571
4572 for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4573 ++i)
4574 temp[i][v] = scratch_data[i][v];
4575 }
4576 }
4577 else
4578 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
4579 template interpolate<false, false>(n_components,
4580 integration_flag,
4581 shape_info,
4582 temp,
4583 values_dofs,
4584 fe_eval.get_face_no());
4585 return false;
4586 }
4587 };
4588
4589
4590
4591 template <int n_face_orientations,
4592 typename Processor,
4593 typename EvaluationData,
4594 const bool check_face_orientations = false>
4595 void
4596 fe_face_evaluation_process_and_io(
4597 Processor & proc,
4598 const unsigned int n_components,
4599 const EvaluationFlags::EvaluationFlags evaluation_flag,
4600 typename Processor::Number2_ * global_vector_ptr,
4601 const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
4602 const EvaluationData & fe_eval,
4603 typename Processor::VectorizedArrayType_ * temp1)
4604 {
4605 constexpr int dim = Processor::dim_;
4606 constexpr int fe_degree = Processor::fe_degree_;
4607 using VectorizedArrayType = typename Processor::VectorizedArrayType_;
4608 constexpr int n_lanes = VectorizedArrayType::size();
4609
4610 using Number = typename Processor::Number_;
4611 using Number2_ = typename Processor::Number2_;
4612
4613 const auto & shape_data = fe_eval.get_shape_info().data.front();
4614 constexpr bool integrate = Processor::do_integrate;
4615 const unsigned int face_no = fe_eval.get_face_no();
4616 const auto & dof_info = fe_eval.get_dof_info();
4617 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
4618 const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
4619 fe_eval.get_dof_access_index();
4620 AssertIndexRange(cell,
4621 dof_info.index_storage_variants[dof_access_index].size());
4622 constexpr unsigned int dofs_per_face =
4623 Utilities::pow(fe_degree + 1, dim - 1);
4624 const unsigned int subface_index = fe_eval.get_subface_index();
4625
4626 const unsigned int n_filled_lanes =
4627 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4628
4629 bool all_faces_are_same = n_filled_lanes == n_lanes;
4630 if (n_face_orientations == n_lanes)
4631 for (unsigned int v = 1; v < n_lanes; ++v)
4632 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
4633 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
4634 {
4635 all_faces_are_same = false;
4636 break;
4637 }
4638
4639 // check for re-orientation ...
4640 std::array<const unsigned int *, n_face_orientations> orientation = {};
4641
4642 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
4643 fe_eval.is_interior_face() == 0)
4644 for (unsigned int v = 0; v < n_lanes; ++v)
4645 {
4646 // the loop breaks once an invalid_unsigned_int is hit for
4647 // all cases except the exterior faces in the ECL loop (where
4648 // some faces might be at the boundaries but others not)
4649 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4650 continue;
4651
4652 if (shape_data.nodal_at_cell_boundaries &&
4653 fe_eval.get_face_orientation(v) != 0)
4654 {
4655 // ... and in case we detect a re-orientation, go to the other
4656 // version of this function that actually allows for this
4657 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4658 check_face_orientations == false)
4659 {
4660 fe_face_evaluation_process_and_io<n_face_orientations,
4661 Processor,
4662 EvaluationData,
4663 true>(proc,
4664 n_components,
4665 evaluation_flag,
4666 global_vector_ptr,
4667 sm_ptr,
4668 fe_eval,
4669 temp1);
4670 return;
4671 }
4672 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4673 fe_eval.get_face_orientation(v), 0);
4674 }
4675 }
4676 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
4677 {
4678 // go to the other version of this function
4679 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4680 check_face_orientations == false)
4681 {
4682 fe_face_evaluation_process_and_io<n_face_orientations,
4683 Processor,
4684 EvaluationData,
4685 true>(proc,
4686 n_components,
4687 evaluation_flag,
4688 global_vector_ptr,
4689 sm_ptr,
4690 fe_eval,
4691 temp1);
4692 return;
4693 }
4694 for (unsigned int v = 0; v < n_face_orientations; ++v)
4695 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4696 fe_eval.get_face_orientation(), 0);
4697 }
4698
4699 // we know that the gradient weights for the Hermite case on the
4700 // right (side==1) are the negative from the value at the left
4701 // (side==0), so we only read out one of them.
4702 VectorizedArrayType grad_weight =
4703 shape_data
4704 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
4705 (1 + face_no % 2))];
4706
4707 // face_to_cell_index_hermite
4708 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
4709 {};
4710 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4711 {
4712 if (n_face_orientations == 1)
4713 index_array_hermite[0] =
4714 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
4715 else
4716 {
4717 for (unsigned int v = 0; v < n_lanes; ++v)
4718 {
4719 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4720 continue;
4721
4722 const auto face_no = fe_eval.get_face_no(v);
4723
4724 grad_weight[v] =
4725 shape_data.shape_data_on_face[0][fe_degree +
4726 (integrate ?
4727 (2 - (face_no % 2)) :
4728 (1 + (face_no % 2)))][0];
4729
4730 index_array_hermite[v] =
4731 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
4732 0);
4733 }
4734 }
4735 }
4736
4737 // face_to_cell_index_nodal
4738 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
4739 {};
4740 if (shape_data.nodal_at_cell_boundaries == true)
4741 {
4742 if (n_face_orientations == 1)
4743 index_array_nodal[0] =
4744 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
4745 else
4746 {
4747 for (unsigned int v = 0; v < n_lanes; ++v)
4748 {
4749 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4750 continue;
4751
4752 const auto face_no = fe_eval.get_face_no(v);
4753
4754 index_array_nodal[v] =
4755 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
4756 0);
4757 }
4758 }
4759 }
4760
4761
4762 const auto reorientate = [&](const unsigned int v, const unsigned int i) {
4763 return (!check_face_orientations || orientation[v] == nullptr) ?
4764 i :
4765 orientation[v][i];
4766 };
4767
4768 const unsigned int cell_index =
4769 dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
4770 fe_eval.get_cell_ids()[0] :
4771 cell * n_lanes;
4772 const unsigned int *dof_indices =
4773 &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
4774
4775 for (unsigned int comp = 0; comp < n_components; ++comp)
4776 {
4777 const std::size_t index_offset =
4778 dof_info.component_dof_indices_offset
4779 [fe_eval.get_active_fe_index()]
4780 [fe_eval.get_first_selected_component()] +
4781 comp * Utilities::pow(fe_degree + 1, dim);
4782
4783 // case 1: contiguous and interleaved indices
4784 if (n_face_orientations == 1 &&
4785 dof_info.index_storage_variants[dof_access_index][cell] ==
4786 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
4787 interleaved_contiguous)
4788 {
4789 AssertDimension(
4790 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4791 n_lanes);
4792 Number2_ *vector_ptr =
4793 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
4794
4795 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4796 {
4797 for (unsigned int i = 0; i < dofs_per_face; ++i)
4798 {
4799 Assert(n_face_orientations == 1, ExcNotImplemented());
4800
4801 const unsigned int ind1 = index_array_hermite[0][2 * i];
4802 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
4803 const unsigned int i_ = reorientate(0, i);
4804 proc.hermite_grad_vectorized(temp1[i_],
4805 temp1[i_ + dofs_per_face],
4806 vector_ptr + ind1 * n_lanes,
4807 vector_ptr + ind2 * n_lanes,
4808 grad_weight);
4809 }
4810 }
4811 else
4812 {
4813 for (unsigned int i = 0; i < dofs_per_face; ++i)
4814 {
4815 Assert(n_face_orientations == 1, ExcNotImplemented());
4816
4817 const unsigned int i_ = reorientate(0, i);
4818 const unsigned int ind = index_array_nodal[0][i];
4819 proc.value_vectorized(temp1[i_],
4820 vector_ptr + ind * n_lanes);
4821 }
4822 }
4823 }
4824
4825 // case 2: contiguous and interleaved indices with fixed stride
4826 else if (n_face_orientations == 1 &&
4827 dof_info.index_storage_variants[dof_access_index][cell] ==
4828 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
4829 interleaved_contiguous_strided)
4830 {
4831 AssertDimension(
4832 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4833 n_lanes);
4834 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
4835 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4836 {
4837 for (unsigned int i = 0; i < dofs_per_face; ++i)
4838 {
4839 Assert(n_face_orientations == 1, ExcNotImplemented());
4840
4841 const unsigned int i_ = reorientate(0, i);
4842 const unsigned int ind1 =
4843 index_array_hermite[0][2 * i] * n_lanes;
4844 const unsigned int ind2 =
4845 index_array_hermite[0][2 * i + 1] * n_lanes;
4846 proc.hermite_grad_vectorized_indexed(
4847 temp1[i_],
4848 temp1[i_ + dofs_per_face],
4849 vector_ptr + ind1,
4850 vector_ptr + ind2,
4851 grad_weight,
4852 dof_indices,
4853 dof_indices);
4854 }
4855 }
4856 else
4857 {
4858 for (unsigned int i = 0; i < dofs_per_face; ++i)
4859 {
4860 Assert(n_face_orientations == 1, ExcNotImplemented());
4861
4862 const unsigned int i_ = reorientate(0, i);
4863 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
4864 proc.value_vectorized_indexed(temp1[i_],
4865 vector_ptr + ind,
4866 dof_indices);
4867 }
4868 }
4869 }
4870
4871 // case 3: contiguous and interleaved indices with mixed stride
4872 else if (n_face_orientations == 1 &&
4873 dof_info.index_storage_variants[dof_access_index][cell] ==
4874 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
4875 interleaved_contiguous_mixed_strides)
4876 {
4877 const unsigned int *strides =
4878 &dof_info.dof_indices_interleave_strides[dof_access_index]
4879 [cell * n_lanes];
4880 unsigned int indices[n_lanes];
4881 for (unsigned int v = 0; v < n_lanes; ++v)
4882 indices[v] = dof_indices[v] + index_offset * strides[v];
4883 const unsigned int n_filled_lanes =
4884 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4885
4886 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4887 {
4888 if (n_filled_lanes == n_lanes)
4889 for (unsigned int i = 0; i < dofs_per_face; ++i)
4890 {
4891 Assert(n_face_orientations == 1, ExcNotImplemented());
4892
4893 const unsigned int i_ = reorientate(0, i);
4894 unsigned int ind1[n_lanes];
4895 DEAL_II_OPENMP_SIMD_PRAGMA
4896 for (unsigned int v = 0; v < n_lanes; ++v)
4897 ind1[v] = indices[v] +
4898 index_array_hermite[0][2 * i] * strides[v];
4899 unsigned int ind2[n_lanes];
4900 DEAL_II_OPENMP_SIMD_PRAGMA
4901 for (unsigned int v = 0; v < n_lanes; ++v)
4902 ind2[v] =
4903 indices[v] +
4904 // TODO
4905 index_array_hermite[0][2 * i + 1] * strides[v];
4906 proc.hermite_grad_vectorized_indexed(
4907 temp1[i_],
4908 temp1[i_ + dofs_per_face],
4909 global_vector_ptr,
4910 global_vector_ptr,
4911 grad_weight,
4912 ind1,
4913 ind2);
4914 }
4915 else
4916 {
4917 if (integrate == false)
4918 for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
4919 temp1[i] = VectorizedArrayType();
4920
4921 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4922 for (unsigned int i = 0; i < dofs_per_face; ++i)
4923 {
4924 const unsigned int i_ =
4925 reorientate(n_face_orientations == 1 ? 0 : v, i);
4926 proc.hermite_grad(
4927 temp1[i_][v],
4928 temp1[i_ + dofs_per_face][v],
4929 global_vector_ptr
4930 [indices[v] +
4931 index_array_hermite
4932 [n_face_orientations == 1 ? 0 : v][2 * i] *
4933 strides[v]],
4934 global_vector_ptr
4935 [indices[v] +
4936 index_array_hermite[n_face_orientations == 1 ?
4937 0 :
4938 v][2 * i + 1] *
4939 strides[v]],
4940 grad_weight[n_face_orientations == 1 ? 0 : v]);
4941 }
4942 }
4943 }
4944 else
4945 {
4946 if (n_filled_lanes == n_lanes)
4947 for (unsigned int i = 0; i < dofs_per_face; ++i)
4948 {
4949 Assert(n_face_orientations == 1, ExcInternalError());
4950 unsigned int ind[n_lanes];
4951 DEAL_II_OPENMP_SIMD_PRAGMA
4952 for (unsigned int v = 0; v < n_lanes; ++v)
4953 ind[v] =
4954 indices[v] + index_array_nodal[0][i] * strides[v];
4955 const unsigned int i_ = reorientate(0, i);
4956 proc.value_vectorized_indexed(temp1[i_],
4957 global_vector_ptr,
4958 ind);
4959 }
4960 else
4961 {
4962 if (integrate == false)
4963 for (unsigned int i = 0; i < dofs_per_face; ++i)
4964 temp1[i] = VectorizedArrayType();
4965
4966 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4967 for (unsigned int i = 0; i < dofs_per_face; ++i)
4968 proc.value(
4969 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
4970 i)][v],
4971 global_vector_ptr
4972 [indices[v] +
4973 index_array_nodal[n_face_orientations == 1 ? 0 : v]
4974 [i] *
4975 strides[v]]);
4976 }
4977 }
4978 }
4979
4980 // case 4: contiguous indices without interleaving
4981 else if (n_face_orientations > 1 ||
4982 dof_info.index_storage_variants[dof_access_index][cell] ==
4983 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
4984 contiguous)
4985 {
4986 Number2_ *vector_ptr = global_vector_ptr + index_offset;
4987
4988 const bool vectorization_possible =
4989 all_faces_are_same && (sm_ptr == nullptr);
4990
4991 std::array<Number2_ *, n_lanes> vector_ptrs;
4992 std::array<unsigned int, n_lanes> reordered_indices;
4993
4994 if (vectorization_possible == false)
4995 {
4996 vector_ptrs = {};
4997 if (n_face_orientations == 1)
4998 {
4999 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5000 if (sm_ptr == nullptr)
5001 {
5002 vector_ptrs[v] = vector_ptr + dof_indices[v];
5003 }
5004 else
5005 {
5006 const auto &temp =
5007 dof_info
5008 .dof_indices_contiguous_sm[dof_access_index]
5009 [cell * n_lanes + v];
5010 vector_ptrs[v] = const_cast<Number2_ *>(
5011 sm_ptr->operator[](temp.first).data() +
5012 temp.second + index_offset);
5013 }
5014 }
5015 else if (n_face_orientations == n_lanes)
5016 {
5017 const auto &cells = fe_eval.get_cell_ids();
5018 for (unsigned int v = 0; v < n_lanes; ++v)
5019 if (cells[v] != numbers::invalid_unsigned_int)
5020 {
5021 if (sm_ptr == nullptr)
5022 {
5023 vector_ptrs[v] =
5024 vector_ptr +
5025 dof_info
5026 .dof_indices_contiguous[dof_access_index]
5027 [cells[v]];
5028 }
5029 else
5030 {
5031 const auto &temp =
5032 dof_info
5033 .dof_indices_contiguous_sm[dof_access_index]
5034 [cells[v]];
5035 vector_ptrs[v] = const_cast<Number2_ *>(
5036 sm_ptr->operator[](temp.first).data() +
5037 temp.second + index_offset);
5038 }
5039 }
5040 }
5041 else
5042 {
5043 Assert(false, ExcNotImplemented());
5044 }
5045 }
5046 else if (n_face_orientations == n_lanes)
5047 {
5048 for (unsigned int v = 0; v < n_lanes; ++v)
5049 reordered_indices[v] =
5050 dof_info.dof_indices_contiguous[dof_access_index]
5051 [fe_eval.get_cell_ids()[v]];
5052 dof_indices = reordered_indices.data();
5053 }
5054
5055 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
5056 {
5057 if (vectorization_possible)
5058 for (unsigned int i = 0; i < dofs_per_face; ++i)
5059 {
5060 const unsigned int ind1 = index_array_hermite[0][2 * i];
5061 const unsigned int ind2 =
5062 index_array_hermite[0][2 * i + 1];
5063 const unsigned int i_ = reorientate(0, i);
5064
5065 proc.hermite_grad_vectorized_indexed(
5066 temp1[i_],
5067 temp1[i_ + dofs_per_face],
5068 vector_ptr + ind1,
5069 vector_ptr + ind2,
5070 grad_weight,
5071 dof_indices,
5072 dof_indices);
5073 }
5074 else if (n_face_orientations == 1)
5075 for (unsigned int i = 0; i < dofs_per_face; ++i)
5076 {
5077 const unsigned int ind1 = index_array_hermite[0][2 * i];
5078 const unsigned int ind2 =
5079 index_array_hermite[0][2 * i + 1];
5080 const unsigned int i_ = reorientate(0, i);
5081
5082 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5083 proc.hermite_grad(temp1[i_][v],
5084 temp1[i_ + dofs_per_face][v],
5085 vector_ptrs[v][ind1],
5086 vector_ptrs[v][ind2],
5087 grad_weight[v]);
5088
5089 if (integrate == false)
5090 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5091 {
5092 temp1[i][v] = 0.0;
5093 temp1[i + dofs_per_face][v] = 0.0;
5094 }
5095 }
5096 else
5097 {
5098 if (integrate == false && n_filled_lanes < n_lanes)
5099 for (unsigned int i = 0; i < dofs_per_face; ++i)
5100 temp1[i] = temp1[i + dofs_per_face] = Number();
5101
5102 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5103 for (unsigned int i = 0; i < dofs_per_face; ++i)
5104 proc.hermite_grad(
5105 temp1[reorientate(v, i)][v],
5106 temp1[reorientate(v, i) + dofs_per_face][v],
5107 vector_ptrs[v][index_array_hermite[v][2 * i]],
5108 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
5109 grad_weight[v]);
5110 }
5111 }
5112 else
5113 {
5114 if (vectorization_possible)
5115 for (unsigned int i = 0; i < dofs_per_face; ++i)
5116 {
5117 const unsigned int ind = index_array_nodal[0][i];
5118 const unsigned int i_ = reorientate(0, i);
5119
5120 proc.value_vectorized_indexed(temp1[i_],
5121 vector_ptr + ind,
5122 dof_indices);
5123 }
5124 else if (n_face_orientations == 1)
5125 for (unsigned int i = 0; i < dofs_per_face; ++i)
5126 {
5127 const unsigned int ind = index_array_nodal[0][i];
5128 const unsigned int i_ = reorientate(0, i);
5129
5130 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5131 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
5132
5133 if (integrate == false)
5134 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5135 temp1[i_][v] = 0.0;
5136 }
5137 else
5138 {
5139 if (integrate == false && n_filled_lanes < n_lanes)
5140 for (unsigned int i = 0; i < dofs_per_face; ++i)
5141 temp1[i] = Number();
5142
5143 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5144 for (unsigned int i = 0; i < dofs_per_face; ++i)
5145 proc.value(temp1[reorientate(v, i)][v],
5146 vector_ptrs[v][index_array_nodal[v][i]]);
5147 }
5148 }
5149 }
5150 else
5151 {
5152 // We should not end up here, this should be caught by
5153 // FEFaceEvaluationImplGatherEvaluateSelector::supports()
5154 Assert(false, ExcInternalError());
5155 }
5156 temp1 += 3 * dofs_per_face;
5157 }
5158 }
5159
5160
5161
5162 template <int dim, typename Number2, typename VectorizedArrayType>
5163 struct FEFaceEvaluationImplGatherEvaluateSelector
5164 {
5165 using Number = typename VectorizedArrayType::value_type;
5166
5167 template <int fe_degree, int n_q_points_1d>
5168 static bool
5169 run(const unsigned int n_components,
5170 const EvaluationFlags::EvaluationFlags evaluation_flag,
5171 const Number2 * src_ptr,
5172 const std::vector<ArrayView<const Number2>> * sm_ptr,
5173 FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
5174 {
5175 Assert(fe_degree > -1, ExcInternalError());
5176 Assert(fe_eval.get_shape_info().element_type <=
5177 MatrixFreeFunctions::tensor_symmetric,
5178 ExcInternalError());
5179
5180 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5181
5182 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5183 VectorizedArrayType *scratch_data =
5184 temp + 3 * n_components * dofs_per_face;
5185
5186 Processor<fe_degree> p;
5187
5188 if (fe_eval.get_dof_access_index() ==
5189 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
5190 fe_eval.is_interior_face() == false)
5191 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5192 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5193 else
5194 fe_face_evaluation_process_and_io<1>(
5195 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5196
5197 const unsigned int subface_index = fe_eval.get_subface_index();
5198
5199 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
5200 FEFaceEvaluationImpl<true,
5201 dim,
5202 fe_degree,
5203 n_q_points_1d,
5204 VectorizedArrayType>::
5205 evaluate_in_face(n_components,
5206 evaluation_flag,
5207 fe_eval.get_shape_info().data.front(),
5208 temp,
5209 fe_eval.begin_values(),
5210 fe_eval.begin_gradients(),
5211 fe_eval.begin_hessians(),
5212 scratch_data,
5213 subface_index);
5214 else
5215 FEFaceEvaluationImpl<false,
5216 dim,
5217 fe_degree,
5218 n_q_points_1d,
5219 VectorizedArrayType>::
5220 evaluate_in_face(n_components,
5221 evaluation_flag,
5222 fe_eval.get_shape_info().data.front(),
5223 temp,
5224 fe_eval.begin_values(),
5225 fe_eval.begin_gradients(),
5226 fe_eval.begin_hessians(),
5227 scratch_data,
5228 subface_index);
5229
5230 // re-orientation for cases not possible with above algorithm
5231 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5232 {
5233 if (fe_eval.get_dof_access_index() ==
5234 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
5235 fe_eval.is_interior_face() == false)
5236 {
5237 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5238 {
5239 // the loop breaks once an invalid_unsigned_int is hit for
5240 // all cases except the exterior faces in the ECL loop (where
5241 // some faces might be at the boundaries but others not)
5242 if (fe_eval.get_cell_ids()[v] ==
5243 numbers::invalid_unsigned_int)
5244 continue;
5245
5246 if (fe_eval.get_face_orientation(v) != 0)
5247 adjust_for_face_orientation_per_lane(
5248 dim,
5249 n_components,
5250 v,
5251 evaluation_flag,
5252 &fe_eval.get_shape_info().face_orientations_quad(
5253 fe_eval.get_face_orientation(v), 0),
5254 false,
5255 Utilities::pow(n_q_points_1d, dim - 1),
5256 &temp[0][0],
5257 fe_eval.begin_values(),
5258 fe_eval.begin_gradients(),
5259 fe_eval.begin_hessians());
5260 }
5261 }
5262 else if (fe_eval.get_face_orientation() != 0)
5263 adjust_for_face_orientation(
5264 dim,
5265 n_components,
5266 evaluation_flag,
5267 &fe_eval.get_shape_info().face_orientations_quad(
5268 fe_eval.get_face_orientation(), 0),
5269 false,
5270 Utilities::pow(n_q_points_1d, dim - 1),
5271 temp,
5272 fe_eval.begin_values(),
5273 fe_eval.begin_gradients(),
5274 fe_eval.begin_hessians());
5275 }
5276
5277 return false;
5278 }
5279
5280 static bool
5281 supports(
5282 const EvaluationFlags::EvaluationFlags evaluation_flag,
5283 const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
5284 const Number2 * vector_ptr,
5285 MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
5286 {
5287 const unsigned int fe_degree = shape_info.data.front().fe_degree;
5288 if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
5289 (evaluation_flag & EvaluationFlags::gradients &&
5290 (fe_degree < 2 ||
5291 shape_info.data.front().element_type !=
5292 MatrixFreeFunctions::tensor_symmetric_hermite)) ||
5293 (evaluation_flag & EvaluationFlags::hessians) ||
5294 vector_ptr == nullptr ||
5295 shape_info.data.front().element_type >
5296 MatrixFreeFunctions::tensor_symmetric ||
5297 storage <
5298 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
5299 return false;
5300 else
5301 return true;
5302 }
5303
5304 private:
5305 template <int fe_degree>
5306 struct Processor
5307 {
5308 static const bool do_integrate = false;
5309 static const int dim_ = dim;
5310 static const int fe_degree_ = fe_degree;
5311 using VectorizedArrayType_ = VectorizedArrayType;
5312 using Number_ = Number;
5313 using Number2_ = const Number2;
5314
5315 template <typename T0, typename T1, typename T2>
5316 void
5317 hermite_grad_vectorized(T0 & temp_1,
5318 T0 & temp_2,
5319 const T1 src_ptr_1,
5320 const T1 src_ptr_2,
5321 const T2 &grad_weight)
5322 {
5323 do_vectorized_read(src_ptr_1, temp_1);
5324 do_vectorized_read(src_ptr_2, temp_2);
5325 temp_2 = grad_weight * (temp_1 - temp_2);
5326 }
5327
5328 template <typename T1, typename T2>
5329 void
5330 value_vectorized(T1 &temp, const T2 src_ptr)
5331 {
5332 do_vectorized_read(src_ptr, temp);
5333 }
5334
5335 template <typename T0, typename T1, typename T2, typename T3>
5336 void
5337 hermite_grad_vectorized_indexed(T0 & temp_1,
5338 T0 & temp_2,
5339 const T1 src_ptr_1,
5340 const T1 src_ptr_2,
5341 const T2 &grad_weight,
5342 const T3 &indices_1,
5343 const T3 &indices_2)
5344 {
5345 do_vectorized_gather(src_ptr_1, indices_1, temp_1);
5346 do_vectorized_gather(src_ptr_2, indices_2, temp_2);
5347 temp_2 = grad_weight * (temp_1 - temp_2);
5348 }
5349
5350 template <typename T0, typename T1, typename T2>
5351 void
5352 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
5353 {
5354 do_vectorized_gather(src_ptr, indices, temp);
5355 }
5356
5357 template <typename T0, typename T1, typename T2>
5358 void
5359 hermite_grad(T0 & temp_1,
5360 T0 & temp_2,
5361 const T1 &src_ptr_1,
5362 const T1 &src_ptr_2,
5363 const T2 &grad_weight)
5364 {
5365 // case 3a)
5366 temp_1 = src_ptr_1;
5367 temp_2 = grad_weight * (temp_1 - src_ptr_2);
5368 }
5369
5370 template <typename T1, typename T2>
5371 void
5372 value(T1 &temp, const T2 &src_ptr)
5373 {
5374 // case 3b)
5375 temp = src_ptr;
5376 }
5377 };
5378 };
5379
5380
5381
5382 template <int dim, typename Number2, typename VectorizedArrayType>
5383 struct FEFaceEvaluationImplIntegrateScatterSelector
5384 {
5385 using Number = typename VectorizedArrayType::value_type;
5386
5387 template <int fe_degree, int n_q_points_1d>
5388 static bool
5389 run(const unsigned int n_components,
5390 const EvaluationFlags::EvaluationFlags integration_flag,
5391 Number2 * dst_ptr,
5392 const std::vector<ArrayView<const Number2>> * sm_ptr,
5393 FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
5394 {
5395 Assert(fe_degree > -1, ExcInternalError());
5396 Assert(fe_eval.get_shape_info().element_type <=
5397 MatrixFreeFunctions::tensor_symmetric,
5398 ExcInternalError());
5399
5400 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5401
5402 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5403 VectorizedArrayType *scratch_data =
5404 temp + 3 * n_components * dofs_per_face;
5405
5406 const unsigned int subface_index = fe_eval.get_subface_index();
5407
5408 // re-orientation for cases not possible with the io function below
5409 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5410 {
5411 if (fe_eval.get_dof_access_index() ==
5412 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
5413 fe_eval.is_interior_face() == false)
5414 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5415 {
5416 // the loop breaks once an invalid_unsigned_int is hit for
5417 // all cases except the exterior faces in the ECL loop (where
5418 // some faces might be at the boundaries but others not)
5419 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
5420 continue;
5421
5422 if (fe_eval.get_face_orientation(v) != 0)
5423 adjust_for_face_orientation_per_lane(
5424 dim,
5425 n_components,
5426 v,
5427 integration_flag,
5428 &fe_eval.get_shape_info().face_orientations_quad(
5429 fe_eval.get_face_orientation(v), 0),
5430 true,
5431 Utilities::pow(n_q_points_1d, dim - 1),
5432 &temp[0][0],
5433 fe_eval.begin_values(),
5434 fe_eval.begin_gradients(),
5435 fe_eval.begin_hessians());
5436 }
5437 else if (fe_eval.get_face_orientation() != 0)
5438 adjust_for_face_orientation(
5439 dim,
5440 n_components,
5441 integration_flag,
5442 &fe_eval.get_shape_info().face_orientations_quad(
5443 fe_eval.get_face_orientation(), 0),
5444 true,
5445 Utilities::pow(n_q_points_1d, dim - 1),
5446 temp,
5447 fe_eval.begin_values(),
5448 fe_eval.begin_gradients(),
5449 fe_eval.begin_hessians());
5450 }
5451
5452 if (fe_degree > -1 && fe_eval.get_subface_index() >=
5453 GeometryInfo<dim - 1>::max_children_per_cell)
5454 FEFaceEvaluationImpl<true,
5455 dim,
5456 fe_degree,
5457 n_q_points_1d,
5458 VectorizedArrayType>::
5459 integrate_in_face(n_components,
5460 integration_flag,
5461 fe_eval.get_shape_info().data.front(),
5462 temp,
5463 fe_eval.begin_values(),
5464 fe_eval.begin_gradients(),
5465 fe_eval.begin_hessians(),
5466 scratch_data,
5467 subface_index);
5468 else
5469 FEFaceEvaluationImpl<false,
5470 dim,
5471 fe_degree,
5472 n_q_points_1d,
5473 VectorizedArrayType>::
5474 integrate_in_face(n_components,
5475 integration_flag,
5476 fe_eval.get_shape_info().data.front(),
5477 temp,
5478 fe_eval.begin_values(),
5479 fe_eval.begin_gradients(),
5480 fe_eval.begin_hessians(),
5481 scratch_data,
5482 subface_index);
5483
5484 Processor<fe_degree> p;
5485
5486 if (fe_eval.get_dof_access_index() ==
5487 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
5488 fe_eval.is_interior_face() == false)
5489 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5490 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5491 else
5492 fe_face_evaluation_process_and_io<1>(
5493 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5494
5495 return false;
5496 }
5497
5498 private:
5499 template <int fe_degree>
5500 struct Processor
5501 {
5502 static const bool do_integrate = true;
5503 static const int dim_ = dim;
5504 static const int fe_degree_ = fe_degree;
5505 using VectorizedArrayType_ = VectorizedArrayType;
5506 using Number_ = Number;
5507 using Number2_ = Number2;
5508
5509 template <typename T0, typename T1, typename T2, typename T3, typename T4>
5510 void
5511 hermite_grad_vectorized(const T0 &temp_1,
5512 const T1 &temp_2,
5513 T2 dst_ptr_1,
5514 T3 dst_ptr_2,
5515 const T4 &grad_weight)
5516 {
5517 // case 1a)
5518 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5519 const VectorizedArrayType grad = grad_weight * temp_2;
5520 do_vectorized_add(val, dst_ptr_1);
5521 do_vectorized_add(grad, dst_ptr_2);
5522 }
5523
5524 template <typename T0, typename T1>
5525 void
5526 value_vectorized(const T0 &temp, T1 dst_ptr)
5527 {
5528 // case 1b)
5529 do_vectorized_add(temp, dst_ptr);
5530 }
5531
5532 template <typename T0, typename T1, typename T2, typename T3>
5533 void
5534 hermite_grad_vectorized_indexed(const T0 &temp_1,
5535 const T0 &temp_2,
5536 T1 dst_ptr_1,
5537 T1 dst_ptr_2,
5538 const T2 &grad_weight,
5539 const T3 &indices_1,
5540 const T3 &indices_2)
5541 {
5542 // case 2a)
5543 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5544 const VectorizedArrayType grad = grad_weight * temp_2;
5545 do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
5546 do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
5547 }
5548
5549 template <typename T0, typename T1, typename T2>
5550 void
5551 value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
5552 {
5553 // case 2b)
5554 do_vectorized_scatter_add(temp, indices, dst_ptr);
5555 }
5556
5557 template <typename T0, typename T1, typename T2>
5558 void
5559 hermite_grad(const T0 &temp_1,
5560 const T0 &temp_2,
5561 T1 & dst_ptr_1,
5562 T1 & dst_ptr_2,
5563 const T2 &grad_weight)
5564 {
5565 // case 3a)
5566 const Number val = temp_1 - grad_weight * temp_2;
5567 const Number grad = grad_weight * temp_2;
5568 dst_ptr_1 += val;
5569 dst_ptr_2 += grad;
5570 }
5571
5572 template <typename T0, typename T1>
5573 void
5574 value(const T0 &temp, T1 &dst_ptr)
5575 {
5576 // case 3b)
5577 dst_ptr += temp;
5578 }
5579 };
5580 };
5581
5582
5583
5588 template <int dim, typename Number>
5589 struct CellwiseInverseMassMatrixImplBasic
5590 {
5591 template <int fe_degree, int = 0>
5592 static bool
5593 run(const unsigned int n_components,
5594 const FEEvaluationData<dim, Number, false> &fe_eval,
5595 const Number * in_array,
5596 Number * out_array,
5597 typename std::enable_if<fe_degree != -1>::type * = nullptr)
5598 {
5599 constexpr unsigned int dofs_per_component =
5600 Utilities::pow(fe_degree + 1, dim);
5601
5602 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5603 Assert(fe_eval.get_shape_info().element_type <=
5604 MatrixFreeFunctions::tensor_symmetric,
5605 ExcNotImplemented());
5606
5607 EvaluatorTensorProduct<evaluate_evenodd,
5608 dim,
5609 fe_degree + 1,
5610 fe_degree + 1,
5611 Number>
5612 evaluator(
5613 AlignedVector<Number>(),
5614 AlignedVector<Number>(),
5615 fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
5616
5617 for (unsigned int d = 0; d < n_components; ++d)
5618 {
5619 const Number *in = in_array + d * dofs_per_component;
5620 Number * out = out_array + d * dofs_per_component;
5621 // Need to select 'apply' method with hessian slot because values
5622 // assume symmetries that do not exist in the inverse shapes
5623 evaluator.template hessians<0, true, false>(in, out);
5624 if (dim > 1)
5625 evaluator.template hessians<1, true, false>(out, out);
5626 if (dim > 2)
5627 evaluator.template hessians<2, true, false>(out, out);
5628 }
5629 for (unsigned int q = 0; q < dofs_per_component; ++q)
5630 {
5631 const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5632 for (unsigned int d = 0; d < n_components; ++d)
5633 out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5634 }
5635 for (unsigned int d = 0; d < n_components; ++d)
5636 {
5637 Number *out = out_array + d * dofs_per_component;
5638 if (dim > 2)
5639 evaluator.template hessians<2, false, false>(out, out);
5640 if (dim > 1)
5641 evaluator.template hessians<1, false, false>(out, out);
5642 evaluator.template hessians<0, false, false>(out, out);
5643 }
5644 return false;
5645 }
5646
5647 template <int fe_degree, int = 0>
5648 static bool
5649 run(const unsigned int n_components,
5650 const FEEvaluationData<dim, Number, false> &fe_eval,
5651 const Number * in_array,
5652 Number * out_array,
5653 typename std::enable_if<fe_degree == -1>::type * = nullptr)
5654 {
5655 static_assert(fe_degree == -1, "Only usable for degree -1");
5656 const unsigned int dofs_per_component =
5657 fe_eval.get_shape_info().dofs_per_component_on_cell;
5658
5659 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5660
5661 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number> evaluator(
5662 fe_eval.get_shape_info().data.front().inverse_shape_values,
5663 AlignedVector<Number>(),
5664 AlignedVector<Number>(),
5665 fe_eval.get_shape_info().data.front().fe_degree + 1,
5666 fe_eval.get_shape_info().data.front().fe_degree + 1);
5667
5668 for (unsigned int d = 0; d < n_components; ++d)
5669 {
5670 const Number *in = in_array + d * dofs_per_component;
5671 Number * out = out_array + d * dofs_per_component;
5672 // Need to select 'apply' method with hessian slot because values
5673 // assume symmetries that do not exist in the inverse shapes
5674 evaluator.template values<0, true, false>(in, out);
5675 if (dim > 1)
5676 evaluator.template values<1, true, false>(out, out);
5677 if (dim > 2)
5678 evaluator.template values<2, true, false>(out, out);
5679 }
5680 for (unsigned int q = 0; q < dofs_per_component; ++q)
5681 {
5682 const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5683 for (unsigned int d = 0; d < n_components; ++d)
5684 out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5685 }
5686 for (unsigned int d = 0; d < n_components; ++d)
5687 {
5688 Number *out = out_array + d * dofs_per_component;
5689 if (dim > 2)
5690 evaluator.template values<2, false, false>(out, out);
5691 if (dim > 1)
5692 evaluator.template values<1, false, false>(out, out);
5693 evaluator.template values<0, false, false>(out, out);
5694 }
5695 return false;
5696 }
5697 };
5698
5699
5700
5705 template <int dim, typename Number>
5706 struct CellwiseInverseMassMatrixImplFlexible
5707 {
5708 template <int fe_degree, int = 0>
5709 static bool
5710 run(const unsigned int n_desired_components,
5711 const AlignedVector<Number> &inverse_shape,
5712 const AlignedVector<Number> &inverse_coefficients,
5713 const Number * in_array,
5714 Number * out_array,
5715 typename std::enable_if<fe_degree != -1>::type * = nullptr)
5716 {
5717 constexpr unsigned int dofs_per_component =
5718 Utilities::pow(fe_degree + 1, dim);
5719 Assert(inverse_coefficients.size() > 0 &&
5720 inverse_coefficients.size() % dofs_per_component == 0,
5721 ExcMessage(
5722 "Expected diagonal to be a multiple of scalar dof per cells"));
5723 if (inverse_coefficients.size() != dofs_per_component)
5724 AssertDimension(n_desired_components * dofs_per_component,
5725 inverse_coefficients.size());
5726
5727 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5728
5729 EvaluatorTensorProduct<evaluate_evenodd,
5730 dim,
5731 fe_degree + 1,
5732 fe_degree + 1,
5733 Number>
5734 evaluator(AlignedVector<Number>(),
5735 AlignedVector<Number>(),
5736 inverse_shape);
5737
5738 const unsigned int shift_coefficient =
5739 inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
5740 0;
5741 const Number *inv_coefficient = inverse_coefficients.data();
5742 for (unsigned int d = 0; d < n_desired_components; ++d)
5743 {
5744 const Number *in = in_array + d * dofs_per_component;
5745 Number * out = out_array + d * dofs_per_component;
5746 // Need to select 'apply' method with hessian slot because values
5747 // assume symmetries that do not exist in the inverse shapes
5748 evaluator.template hessians<0, true, false>(in, out);
5749 if (dim > 1)
5750 evaluator.template hessians<1, true, false>(out, out);
5751 if (dim > 2)
5752 evaluator.template hessians<2, true, false>(out, out);
5753
5754 for (unsigned int q = 0; q < dofs_per_component; ++q)
5755 out[q] *= inv_coefficient[q];
5756
5757 if (dim > 2)
5758 evaluator.template hessians<2, false, false>(out, out);
5759 if (dim > 1)
5760 evaluator.template hessians<1, false, false>(out, out);
5761 evaluator.template hessians<0, false, false>(out, out);
5762
5763 inv_coefficient += shift_coefficient;
5764 }
5765 return false;
5766 }
5767
5771 template <int fe_degree, int = 0>
5772 static bool
5773 run(const unsigned int,
5774 const AlignedVector<Number> &,
5775 const AlignedVector<Number> &,
5776 const Number *,
5777 Number *,
5778 typename std::enable_if<fe_degree == -1>::type * = nullptr)
5779 {
5780 static_assert(fe_degree == -1, "Only usable for degree -1");
5781 Assert(false, ExcNotImplemented());
5782 return false;
5783 }
5784 };
5785
5786
5787
5792 template <int dim, typename Number>
5793 struct CellwiseInverseMassMatrixImplTransformFromQPoints
5794 {
5795 template <int fe_degree, int n_q_points_1d>
5796 static bool
5797 run(const unsigned int n_desired_components,
5798 const FEEvaluationData<dim, Number, false> &fe_eval,
5799 const Number * in_array,
5800 Number * out_array)
5801 {
5802 static const bool do_inplace =
5803 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
5804
5805 Assert(fe_eval.get_shape_info().element_type !=
5806 MatrixFreeFunctions::ElementType::tensor_none,
5807 ExcNotImplemented());
5808
5809 const auto &inverse_shape =
5810 do_inplace ?
5811 fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
5812 fe_eval.get_shape_info().data.front().inverse_shape_values;
5813
5814 const std::size_t dofs_per_component =
5815 do_inplace ? Utilities::pow(fe_degree + 1, dim) :
5816 fe_eval.get_shape_info().dofs_per_component_on_cell;
5817 const std::size_t n_q_points = do_inplace ?
5818 Utilities::pow(fe_degree + 1, dim) :
5819 fe_eval.get_shape_info().n_q_points;
5820
5821 EvaluatorTensorProduct<do_inplace ? evaluate_evenodd : evaluate_general,
5822 dim,
5823 fe_degree + 1,
5824 n_q_points_1d,
5825 Number>
5826 evaluator(AlignedVector<Number>(),
5827 AlignedVector<Number>(),
5828 inverse_shape,
5829 fe_eval.get_shape_info().data.front().fe_degree + 1,
5830 fe_eval.get_shape_info().data.front().n_q_points_1d);
5831
5832 for (unsigned int d = 0; d < n_desired_components; ++d)
5833 {
5834 const Number *in = in_array + d * n_q_points;
5835 Number * out = out_array + d * dofs_per_component;
5836
5837 auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
5838 auto temp_2 = do_inplace ?
5839 out :
5840 (temp_1 + std::max(n_q_points, dofs_per_component));
5841
5842 if (dim == 3)
5843 {
5844 evaluator.template hessians<2, false, false>(in, temp_1);
5845 evaluator.template hessians<1, false, false>(temp_1, temp_2);
5846 evaluator.template hessians<0, false, false>(temp_2, out);
5847 }
5848 if (dim == 2)
5849 {
5850 evaluator.template hessians<1, false, false>(in, temp_1);
5851 evaluator.template hessians<0, false, false>(temp_1, out);
5852 }
5853 if (dim == 1)
5854 evaluator.template hessians<0, false, false>(in, out);
5855 }
5856 return false;
5857 }
5858 };
5859
5860} // end of namespace internal
5861
5862
5863DEAL_II_NAMESPACE_CLOSE
5864
5865#endif
AlignedVector::data
pointer data()
AlignedVector::size
size_type size() const
ArrayView::begin
iterator begin() const
Definition: array_view.h:585
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::get_shape_info
const ShapeInfoType & get_shape_info() const
FEEvaluationData::JxW
Number JxW(const unsigned int q_point) const
FEEvaluationData::get_cell_ids
const std::array< unsigned int, n_lanes > & get_cell_ids() const
Definition: fe_evaluation_data.h:481
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:628
VectorizedArray::load
void load(const Number *ptr)
Definition: vectorization.h:543
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
evaluation_flags.h
fe_evaluation_data.h
geometry_info.h
cell_index
unsigned int cell_index
Definition: grid_tools.cc:1129
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
internal::MatrixFreeFunctions::ElementType
ElementType
Definition: shape_info.h:59
internal::MatrixFreeFunctions::tensor_raviart_thomas
@ tensor_raviart_thomas
Definition: shape_info.h:107
internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0
@ tensor_symmetric_plus_dg0
Definition: shape_info.h:100
internal::MatrixFreeFunctions::truncated_tensor
@ truncated_tensor
Definition: shape_info.h:93
internal::MatrixFreeFunctions::tensor_none
@ tensor_none
Definition: shape_info.h:112
internal::MatrixFreeFunctions::tensor_symmetric
@ tensor_symmetric
Definition: shape_info.h:82
internal::MatrixFreeFunctions::tensor_symmetric_hermite
@ tensor_symmetric_hermite
Definition: shape_info.h:75
dof_info.h
EvaluationFlags::EvaluationFlags
EvaluationFlags
The EvaluationFlags enum.
Definition: evaluation_flags.h:43
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1123
internal::do_vectorized_add
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
Definition: evaluation_kernels.h:3927
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition: evaluation_kernels.h:2243
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.h:4046
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.h:4596
internal::do_vectorized_scatter_add
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
Definition: evaluation_kernels.h:3951
internal::do_vectorized_gather
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
Definition: evaluation_kernels.h:3901
internal::do_vectorized_read
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
Definition: evaluation_kernels.h:3879
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.h:3983
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5605
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5588
shape_info.h
internal::CellwiseInverseMassMatrixImplBasic::run
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
Definition: evaluation_kernels.h:5593
internal::CellwiseInverseMassMatrixImplBasic::run
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
Definition: evaluation_kernels.h:5649
internal::CellwiseInverseMassMatrixImplFlexible::run
static bool run(const unsigned int n_desired_components, const AlignedVector< Number > &inverse_shape, const AlignedVector< Number > &inverse_coefficients, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
Definition: evaluation_kernels.h:5710
internal::CellwiseInverseMassMatrixImplFlexible::run
static bool run(const unsigned int, const AlignedVector< Number > &, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
Definition: evaluation_kernels.h:5773
internal::CellwiseInverseMassMatrixImplTransformFromQPoints::run
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
Definition: evaluation_kernels.h:5797
internal::FEEvaluationImplBasisChange::do_mass
static void do_mass(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
Definition: evaluation_kernels.h:1742
internal::FEEvaluationImplBasisChange::do_backward
static void do_backward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const bool add_into_result, Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
Definition: evaluation_kernels.h:1592
internal::FEEvaluationImplBasisChange::do_forward
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
Definition: evaluation_kernels.h:1468
internal::FEEvaluationImplCollocation::integrate
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
Definition: evaluation_kernels.h:1958
internal::FEEvaluationImplCollocation::do_integrate
static void do_integrate(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, Number *gradients_quad, const Number *hessians_quad, const bool add_into_values_array)
Definition: evaluation_kernels.h:1996
internal::FEEvaluationImplCollocation::do_evaluate
static void do_evaluate(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, Number *gradients_quad, Number *hessians_quad)
Definition: evaluation_kernels.h:1897
internal::FEEvaluationImplCollocation::evaluate
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition: evaluation_kernels.h:1869
internal::FEEvaluationImplEvaluateSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition: evaluation_kernels.h:2271
internal::FEEvaluationImplIntegrateSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
Definition: evaluation_kernels.h:2400
internal::FEEvaluationImplTransformToCollocation::integrate
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
Definition: evaluation_kernels.h:2185
internal::FEEvaluationImplTransformToCollocation::evaluate
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition: evaluation_kernels.h:2134
internal::FEEvaluationImpl< MatrixFreeFunctions::tensor_raviart_thomas, dim, fe_degree, n_q_points_1d, Number >::create_evaluator_tensor_product
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape_data)
Definition: evaluation_kernels.h:210
internal::FEEvaluationImpl::variant
static const EvaluatorVariant variant
Definition: evaluation_kernels.h:119
internal::FEEvaluationImpl::create_evaluator_tensor_product
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > *univariate_shape_data)
Definition: evaluation_kernels.h:142
internal::FEEvaluationImpl::integrate
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
Definition: evaluation_kernels.h:518
internal::FEEvaluationImpl::Eval
EvaluatorTensorProduct< variant, dim, fe_degree+1, n_q_points_1d, Number > Eval
Definition: evaluation_kernels.h:126
internal::FEEvaluationImpl::evaluate
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval)
Definition: evaluation_kernels.h:245
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.h:4117
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.h:5359
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::value_vectorized_indexed
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
Definition: evaluation_kernels.h:5352
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.h:5317
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.h:5337
internal::FEFaceEvaluationImplGatherEvaluateSelector::Number
typename VectorizedArrayType::value_type Number
Definition: evaluation_kernels.h:5165
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.h:5169
internal::FEFaceEvaluationImplGatherEvaluateSelector::supports
static bool supports(const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, const Number2 *vector_ptr, MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
Definition: evaluation_kernels.h:5281
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.h:5511
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.h:5534
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.h:5559
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::value_vectorized_indexed
void value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
Definition: evaluation_kernels.h:5551
internal::FEFaceEvaluationImplIntegrateScatterSelector::Number
typename VectorizedArrayType::value_type Number
Definition: evaluation_kernels.h:5385
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.h:5389
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)
Definition: evaluation_kernels.h:4358
internal::FEFaceEvaluationImplRaviartThomas::evaluate_or_integrate_in_face
static void evaluate_or_integrate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, Number *scratch_data, const unsigned int subface_index, const unsigned int face_no)
Definition: evaluation_kernels.h:3037
internal::FEFaceEvaluationImplRaviartThomas::evaluate_in_face_apply
static void evaluate_in_face_apply(Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, Number *scratch_data, const EvaluationFlags::EvaluationFlags evaluation_flag, const unsigned int face_direction, const unsigned int subface_index, std::integral_constant< bool, false >)
Definition: evaluation_kernels.h:3101
internal::FEFaceEvaluationImplRaviartThomas::evaluate_in_face_apply
static void evaluate_in_face_apply(Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, Number *scratch_data, const EvaluationFlags::EvaluationFlags evaluation_flag, const unsigned int face_direction, const unsigned int subface_index, std::integral_constant< bool, true >)
Definition: evaluation_kernels.h:3305
internal::FEFaceEvaluationImplRaviartThomas::create_evaluator_tensor_product
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > &data, const unsigned int subface_index, const unsigned int direction)
Definition: evaluation_kernels.h:3016
internal::FEFaceEvaluationImpl::evaluate_in_face
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::UnivariateShapeData< Number > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
Definition: evaluation_kernels.h:2571
internal::FEFaceEvaluationImpl::create_evaluator_tensor_product
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > &data, const unsigned int subface_index, const unsigned int direction)
Definition: evaluation_kernels.h:2541
internal::FEFaceEvaluationImpl::integrate_in_face
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::UnivariateShapeData< Number > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
Definition: evaluation_kernels.h:2772
internal::FEFaceEvaluationImpl::Eval
EvaluatorTensorProduct< symmetric_evaluate ? evaluate_evenodd :evaluate_general, dim - 1, fe_degree+1, n_q_points_1d, Number > Eval
Definition: evaluation_kernels.h:2538
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< Number >, 2 > &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face)
Definition: evaluation_kernels.h:3604
internal::FEFaceNormalEvaluationImpl::interpolate
static void interpolate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *input, Number *output, const unsigned int face_no)
Definition: evaluation_kernels.h:3544
internal::FEFaceNormalEvaluationImpl::interpolate_generic_raviart_thomas_apply_face
static void interpolate_generic_raviart_thomas_apply_face(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const unsigned int face_no, const Number *input, Number *output)
Definition: evaluation_kernels.h:3756
internal::FEFaceNormalEvaluationImpl::create_evaluator_tensor_product
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > &data, const unsigned int face_no)
Definition: evaluation_kernels.h:3674
internal::FEFaceNormalEvaluationImpl::interpolate_generic_raviart_thomas
static void interpolate_generic_raviart_thomas(const unsigned int n_components, const Number *input, Number *output, const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info)
Definition: evaluation_kernels.h:3688
internal::FEFaceNormalEvaluationImpl::interpolate_quadrature
static void interpolate_quadrature(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *input, Number *output, const unsigned int face_no)
Definition: evaluation_kernels.h:3576
internal::MatrixFreeFunctions::ShapeInfo::data
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:423
internal::MatrixFreeFunctions::ShapeInfo::face_orientations_quad
::Table< 2, unsigned int > face_orientations_quad
Definition: shape_info.h:562
internal::MatrixFreeFunctions::UnivariateShapeData::shape_data_on_face
std::array< AlignedVector< Number >, 2 > shape_data_on_face
Definition: shape_info.h:247
internal::MatrixFreeFunctions::UnivariateShapeData::shape_gradients_collocation_eo
AlignedVector< Number > shape_gradients_collocation_eo
Definition: shape_info.h:211
internal::MatrixFreeFunctions::UnivariateShapeData::hessians_within_subface
std::array< AlignedVector< Number >, 2 > hessians_within_subface
Definition: shape_info.h:277
internal::MatrixFreeFunctions::UnivariateShapeData::values_within_subface
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition: shape_info.h:265
internal::MatrixFreeFunctions::UnivariateShapeData::shape_hessians_collocation_eo
AlignedVector< Number > shape_hessians_collocation_eo
Definition: shape_info.h:219
internal::MatrixFreeFunctions::UnivariateShapeData::gradients_within_subface
std::array< AlignedVector< Number >, 2 > gradients_within_subface
Definition: shape_info.h:271
tensor_product_kernels.h
Generated by doxygen 1.9.5