deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
evaluation_kernels.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2017 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_evaluation_kernels_h
15#define dealii_matrix_free_evaluation_kernels_h
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/utilities.h>
20#include <deal.II/base/vectorization.h>
21
22#include <deal.II/matrix_free/evaluation_flags.h>
23#include <deal.II/matrix_free/evaluation_kernels_common.h>
24#include <deal.II/matrix_free/fe_evaluation_data.h>
25#include <deal.II/matrix_free/shape_info.h>
26#include <deal.II/matrix_free/tensor_product_kernels.h>
27
28
29DEAL_II_NAMESPACE_OPEN
30
31
32namespace internal
33{
34 // Select evaluator type from element shape function type
35 template <MatrixFreeFunctions::ElementType element, bool is_long>
36 struct EvaluatorSelector
37 {};
38
39 template <bool is_long>
40 struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
41 {
42 static const EvaluatorVariant variant = evaluate_general;
43 };
44
45 template <>
46 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
47 {
48 static const EvaluatorVariant variant = evaluate_symmetric;
49 };
50
51 template <>
52 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
53 {
54 static const EvaluatorVariant variant = evaluate_evenodd;
55 };
56
57 template <bool is_long>
58 struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
59 {
60 static const EvaluatorVariant variant = evaluate_general;
61 };
62
63 template <>
64 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
65 false>
66 {
67 static const EvaluatorVariant variant = evaluate_general;
68 };
69
70 template <>
71 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
72 {
73 static const EvaluatorVariant variant = evaluate_evenodd;
74 };
75
76 template <bool is_long>
77 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
78 is_long>
79 {
80 static const EvaluatorVariant variant = evaluate_evenodd;
81 };
82
83
84
106 template <MatrixFreeFunctions::ElementType type,
107 int dim,
108 int fe_degree,
109 int n_q_points_1d,
110 typename Number>
111 struct FEEvaluationImpl
112 {
113 static const EvaluatorVariant variant =
114 EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
115 using Number2 =
116 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
117
118 using Eval = EvaluatorTensorProduct<variant,
119 dim,
120 fe_degree + 1,
121 n_q_points_1d,
122 Number,
123 Number2>;
124
125 static void
126 evaluate(const unsigned int n_components,
127 const EvaluationFlags::EvaluationFlags evaluation_flag,
128 const Number *values_dofs_actual,
129 FEEvaluationData<dim, Number, false> &fe_eval);
130
131 static void
132 integrate(const unsigned int n_components,
133 const EvaluationFlags::EvaluationFlags integration_flag,
134 Number *values_dofs_actual,
135 FEEvaluationData<dim, Number, false> &fe_eval,
136 const bool add_into_values_array);
137
138 static Eval
139 create_evaluator_tensor_product(
140 const MatrixFreeFunctions::UnivariateShapeData<Number2>
141 *univariate_shape_data)
142 {
143 if (variant == evaluate_evenodd)
144 return Eval(univariate_shape_data->shape_values_eo,
145 univariate_shape_data->shape_gradients_eo,
146 univariate_shape_data->shape_hessians_eo,
147 univariate_shape_data->fe_degree + 1,
148 univariate_shape_data->n_q_points_1d);
149 else
150 return Eval(univariate_shape_data->shape_values,
151 univariate_shape_data->shape_gradients,
152 univariate_shape_data->shape_hessians,
153 univariate_shape_data->fe_degree + 1,
154 univariate_shape_data->n_q_points_1d);
155 }
156 };
157
158
159
164 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
165 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
166 dim,
167 fe_degree,
168 n_q_points_1d,
169 Number>
170 {
171 static void
172 evaluate(const unsigned int n_components,
173 const EvaluationFlags::EvaluationFlags evaluation_flag,
174 const Number *values_dofs_actual,
175 FEEvaluationData<dim, Number, false> &fe_eval);
176
177 static void
178 integrate(const unsigned int n_components,
179 const EvaluationFlags::EvaluationFlags integration_flag,
180 Number *values_dofs_actual,
181 FEEvaluationData<dim, Number, false> &fe_eval,
182 const bool add_into_values_array);
183 };
184
185
186
187 template <MatrixFreeFunctions::ElementType type,
188 int dim,
189 int fe_degree,
190 int n_q_points_1d,
191 typename Number>
192 inline void
193 FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::evaluate(
194 const unsigned int n_components,
195 const EvaluationFlags::EvaluationFlags evaluation_flag,
196 const Number *values_dofs_actual,
197 FEEvaluationData<dim, Number, false> &fe_eval)
198 {
199 if (evaluation_flag == EvaluationFlags::nothing)
200 return;
201
202 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
203 univariate_shape_data;
204
205 const auto &shape_data = fe_eval.get_shape_info().data;
206
207 univariate_shape_data.fill(&shape_data.front());
208
209 if (shape_data.size() == dim)
210 for (int i = 1; i < dim; ++i)
211 univariate_shape_data[i] = &shape_data[i];
212
213 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
214 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
215 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
216
217 const unsigned int temp_size =
218 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
219 0 :
220 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
221 Eval::n_rows_of_product :
222 Eval::n_columns_of_product);
223 Number *temp1 = fe_eval.get_scratch_data().begin();
224 Number *temp2;
225 if (temp_size == 0)
226 {
227 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
228 shape_data.front().fe_degree + 1),
229 Utilities::fixed_power<dim>(
230 shape_data.front().n_q_points_1d));
231 }
232 else
233 {
234 temp2 = temp1 + temp_size;
235 }
236
237 const std::size_t n_q_points = temp_size == 0 ?
238 fe_eval.get_shape_info().n_q_points :
239 Eval::n_columns_of_product;
240 const std::size_t dofs_per_comp =
241 (type == MatrixFreeFunctions::truncated_tensor) ?
242 Utilities::pow(shape_data.front().fe_degree + 1, dim) :
243 fe_eval.get_shape_info().dofs_per_component_on_cell;
244 const Number *values_dofs =
245 (type == MatrixFreeFunctions::truncated_tensor) ?
246 temp1 + 2 * (std::max<std::size_t>(
247 fe_eval.get_shape_info().dofs_per_component_on_cell,
248 n_q_points)) :
249 values_dofs_actual;
250
251 if (type == MatrixFreeFunctions::truncated_tensor)
252 embed_truncated_into_full_tensor_product<dim, fe_degree>(
253 n_components,
254 const_cast<Number *>(values_dofs),
255 values_dofs_actual,
256 fe_eval);
257
258 Number *values_quad = fe_eval.begin_values();
259 Number *gradients_quad = fe_eval.begin_gradients();
260
261 switch (dim)
262 {
263 case 1:
264 for (unsigned int c = 0; c < n_components; ++c)
265 {
266 if (evaluation_flag & EvaluationFlags::values)
267 eval0.template values<0, true, false>(values_dofs, values_quad);
268 if (evaluation_flag & EvaluationFlags::gradients)
269 eval0.template gradients<0, true, false>(values_dofs,
270 gradients_quad);
271
272 // advance the next component in 1d array
273 values_dofs += dofs_per_comp;
274 values_quad += n_q_points;
275 gradients_quad += n_q_points;
276 }
277 break;
278
279 case 2:
280 for (unsigned int c = 0; c < n_components; ++c)
281 {
282 // grad x
283 if (evaluation_flag & EvaluationFlags::gradients)
284 {
285 eval0.template gradients<0, true, false>(values_dofs, temp1);
286 eval1.template values<1, true, false, 2>(temp1,
287 gradients_quad);
288 }
289
290 // grad y
291 eval0.template values<0, true, false>(values_dofs, temp1);
292 if (evaluation_flag & EvaluationFlags::gradients)
293 eval1.template gradients<1, true, false, 2>(temp1,
294 gradients_quad + 1);
295
296 // val: can use values applied in x
297 if (evaluation_flag & EvaluationFlags::values)
298 eval1.template values<1, true, false>(temp1, values_quad);
299
300 // advance to the next component in 1d array
301 values_dofs += dofs_per_comp;
302 values_quad += n_q_points;
303 gradients_quad += 2 * n_q_points;
304 }
305 break;
306
307 case 3:
308 for (unsigned int c = 0; c < n_components; ++c)
309 {
310 if (evaluation_flag & EvaluationFlags::gradients)
311 {
312 // grad x
313 eval0.template gradients<0, true, false>(values_dofs, temp1);
314 eval1.template values<1, true, false>(temp1, temp2);
315 eval2.template values<2, true, false, 3>(temp2,
316 gradients_quad);
317 }
318
319 // grad y
320 eval0.template values<0, true, false>(values_dofs, temp1);
321 if (evaluation_flag & EvaluationFlags::gradients)
322 {
323 eval1.template gradients<1, true, false>(temp1, temp2);
324 eval2.template values<2, true, false, 3>(temp2,
325 gradients_quad + 1);
326 }
327
328 // grad z: can use the values applied in x direction stored in
329 // temp1
330 eval1.template values<1, true, false>(temp1, temp2);
331 if (evaluation_flag & EvaluationFlags::gradients)
332 eval2.template gradients<2, true, false, 3>(temp2,
333 gradients_quad + 2);
334
335 // val: can use the values applied in x & y direction stored in
336 // temp2
337 if (evaluation_flag & EvaluationFlags::values)
338 eval2.template values<2, true, false>(temp2, values_quad);
339
340 // advance to the next component in 1d array
341 values_dofs += dofs_per_comp;
342 values_quad += n_q_points;
343 gradients_quad += 3 * n_q_points;
344 }
345 break;
346
347 default:
348 AssertThrow(false, ExcNotImplemented());
349 }
350
351 // case additional dof for FE_Q_DG0: add values; gradients and second
352 // derivatives evaluate to zero
353 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 &&
354 (evaluation_flag & EvaluationFlags::values))
355 {
356 values_quad -= n_components * n_q_points;
357 values_dofs -= n_components * dofs_per_comp;
358 for (std::size_t c = 0; c < n_components; ++c)
359 for (std::size_t q = 0; q < n_q_points; ++q)
360 values_quad[c * n_q_points + q] +=
361 values_dofs[(c + 1) * dofs_per_comp - 1];
362 }
363 }
364
365
366
367 template <MatrixFreeFunctions::ElementType type,
368 int dim,
369 int fe_degree,
370 int n_q_points_1d,
371 typename Number>
372 inline void
373 FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::integrate(
374 const unsigned int n_components,
375 const EvaluationFlags::EvaluationFlags integration_flag,
376 Number *values_dofs_actual,
377 FEEvaluationData<dim, Number, false> &fe_eval,
378 const bool add_into_values_array)
379 {
380 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
381 univariate_shape_data;
382
383 const auto &shape_data = fe_eval.get_shape_info().data;
384 univariate_shape_data.fill(&shape_data.front());
385
386 if (shape_data.size() == dim)
387 for (int i = 1; i < dim; ++i)
388 univariate_shape_data[i] = &shape_data[i];
389
390 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
391 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
392 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
393
394 const unsigned int temp_size =
395 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
396 0 :
397 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
398 Eval::n_rows_of_product :
399 Eval::n_columns_of_product);
400 Number *temp1 = fe_eval.get_scratch_data().begin();
401 Number *temp2;
402 if (temp_size == 0)
403 {
404 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
405 shape_data.front().fe_degree + 1),
406 Utilities::fixed_power<dim>(
407 shape_data.front().n_q_points_1d));
408 }
409 else
410 {
411 temp2 = temp1 + temp_size;
412 }
413
414 const std::size_t n_q_points = temp_size == 0 ?
415 fe_eval.get_shape_info().n_q_points :
416 Eval::n_columns_of_product;
417 const unsigned int dofs_per_comp =
418 (type == MatrixFreeFunctions::truncated_tensor) ?
419 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
420 fe_eval.get_shape_info().dofs_per_component_on_cell;
421
422 // expand dof_values to tensor product for truncated tensor products
423 Number *values_dofs =
424 (type == MatrixFreeFunctions::truncated_tensor) ?
425 temp1 + 2 * (std::max<std::size_t>(
426 fe_eval.get_shape_info().dofs_per_component_on_cell,
427 n_q_points)) :
428 values_dofs_actual;
429
430 Number *values_quad = fe_eval.begin_values();
431 Number *gradients_quad = fe_eval.begin_gradients();
432
433 switch (dim)
434 {
435 case 1:
436 for (unsigned int c = 0; c < n_components; ++c)
437 {
438 if (integration_flag & EvaluationFlags::values)
439 {
440 if (add_into_values_array == false)
441 eval0.template values<0, false, false>(values_quad,
442 values_dofs);
443 else
444 eval0.template values<0, false, true>(values_quad,
445 values_dofs);
446 }
447 if (integration_flag & EvaluationFlags::gradients)
448 {
449 if (integration_flag & EvaluationFlags::values ||
450 add_into_values_array == true)
451 eval0.template gradients<0, false, true>(gradients_quad,
452 values_dofs);
453 else
454 eval0.template gradients<0, false, false>(gradients_quad,
455 values_dofs);
456 }
457
458 // advance to the next component in 1d array
459 values_dofs += dofs_per_comp;
460 values_quad += n_q_points;
461 gradients_quad += n_q_points;
462 }
463 break;
464
465 case 2:
466 for (unsigned int c = 0; c < n_components; ++c)
467 {
468 if ((integration_flag & EvaluationFlags::values) &&
469 !(integration_flag & EvaluationFlags::gradients))
470 {
471 eval1.template values<1, false, false>(values_quad, temp1);
472 if (add_into_values_array == false)
473 eval0.template values<0, false, false>(temp1, values_dofs);
474 else
475 eval0.template values<0, false, true>(temp1, values_dofs);
476 }
477 if (integration_flag & EvaluationFlags::gradients)
478 {
479 eval1.template gradients<1, false, false, 2>(gradients_quad +
480 1,
481 temp1);
482 if (integration_flag & EvaluationFlags::values)
483 eval1.template values<1, false, true>(values_quad, temp1);
484 if (add_into_values_array == false)
485 eval0.template values<0, false, false>(temp1, values_dofs);
486 else
487 eval0.template values<0, false, true>(temp1, values_dofs);
488 eval1.template values<1, false, false, 2>(gradients_quad,
489 temp1);
490 eval0.template gradients<0, false, true>(temp1, values_dofs);
491 }
492
493 // advance to the next component in 1d array
494 values_dofs += dofs_per_comp;
495 values_quad += n_q_points;
496 gradients_quad += 2 * n_q_points;
497 }
498 break;
499
500 case 3:
501 for (unsigned int c = 0; c < n_components; ++c)
502 {
503 if ((integration_flag & EvaluationFlags::values) &&
504 !(integration_flag & EvaluationFlags::gradients))
505 {
506 eval2.template values<2, false, false>(values_quad, temp1);
507 eval1.template values<1, false, false>(temp1, temp2);
508 if (add_into_values_array == false)
509 eval0.template values<0, false, false>(temp2, values_dofs);
510 else
511 eval0.template values<0, false, true>(temp2, values_dofs);
512 }
513 if (integration_flag & EvaluationFlags::gradients)
514 {
515 eval2.template gradients<2, false, false, 3>(gradients_quad +
516 2,
517 temp1);
518 if (integration_flag & EvaluationFlags::values)
519 eval2.template values<2, false, true>(values_quad, temp1);
520 eval1.template values<1, false, false>(temp1, temp2);
521 eval2.template values<2, false, false, 3>(gradients_quad + 1,
522 temp1);
523 eval1.template gradients<1, false, true>(temp1, temp2);
524 if (add_into_values_array == false)
525 eval0.template values<0, false, false>(temp2, values_dofs);
526 else
527 eval0.template values<0, false, true>(temp2, values_dofs);
528 eval2.template values<2, false, false, 3>(gradients_quad,
529 temp1);
530 eval1.template values<1, false, false>(temp1, temp2);
531 eval0.template gradients<0, false, true>(temp2, values_dofs);
532 }
533
534 // advance to the next component in 1d array
535 values_dofs += dofs_per_comp;
536 values_quad += n_q_points;
537 gradients_quad += 3 * n_q_points;
538 }
539 break;
540
541 default:
542 AssertThrow(false, ExcNotImplemented());
543 }
544
545 // case FE_Q_DG0: add values, gradients and second derivatives are zero
546 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
547 {
548 values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
549 values_quad -= n_components * n_q_points;
550 if (integration_flag & EvaluationFlags::values)
551 for (unsigned int c = 0; c < n_components; ++c)
552 {
553 values_dofs[0] = values_quad[0];
554 for (unsigned int q = 1; q < n_q_points; ++q)
555 values_dofs[0] += values_quad[q];
556 values_dofs += dofs_per_comp;
557 values_quad += n_q_points;
558 }
559 else
560 {
561 for (unsigned int c = 0; c < n_components; ++c)
562 values_dofs[c * dofs_per_comp] = Number();
563 values_dofs += n_components * dofs_per_comp;
564 }
565 }
566
567 if (type == MatrixFreeFunctions::truncated_tensor)
568 truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
569 n_components,
570 values_dofs_actual,
571 values_dofs - dofs_per_comp * n_components,
572 fe_eval);
573 }
574
575
576
577 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
578 inline void
579 FEEvaluationImpl<
580 MatrixFreeFunctions::tensor_none,
581 dim,
582 fe_degree,
583 n_q_points_1d,
584 Number>::evaluate(const unsigned int n_components,
585 const EvaluationFlags::EvaluationFlags evaluation_flag,
586 const Number *values_dofs_actual,
587 FEEvaluationData<dim, Number, false> &fe_eval)
588 {
589 Assert(!(evaluation_flag & EvaluationFlags::hessians), ExcNotImplemented());
590
591 const std::size_t n_dofs =
592 fe_eval.get_shape_info().dofs_per_component_on_cell;
593 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
594
595 const auto &shape_data = fe_eval.get_shape_info().data;
596
597 using Number2 =
598 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
599
600 if (evaluation_flag & EvaluationFlags::values)
601 {
602 const auto *const shape_values = shape_data.front().shape_values.data();
603 auto *out = fe_eval.begin_values();
604 const auto *in = values_dofs_actual;
605
606 for (unsigned int c = 0; c < n_components; c += 3)
607 {
608 if (c + 1 == n_components)
609 apply_matrix_vector_product<evaluate_general,
610 EvaluatorQuantity::value,
611 /*transpose_matrix*/ true,
612 /*add*/ false,
613 /*consider_strides*/ false,
614 Number,
615 Number2,
616 /*n_components*/ 1>(
617 shape_values, in, out, n_dofs, n_q_points, 1, 1);
618 else if (c + 2 == n_components)
619 apply_matrix_vector_product<evaluate_general,
620 EvaluatorQuantity::value,
621 /*transpose_matrix*/ true,
622 /*add*/ false,
623 /*consider_strides*/ false,
624 Number,
625 Number2,
626 /*n_components*/ 2>(
627 shape_values, in, out, n_dofs, n_q_points, 1, 1);
628 else
629 apply_matrix_vector_product<evaluate_general,
630 EvaluatorQuantity::value,
631 /*transpose_matrix*/ true,
632 /*add*/ false,
633 /*consider_strides*/ false,
634 Number,
635 Number2,
636 /*n_components*/ 3>(
637 shape_values, in, out, n_dofs, n_q_points, 1, 1);
638
639 out += 3 * n_q_points;
640 in += 3 * n_dofs;
641 }
642 }
643
644 if (evaluation_flag & EvaluationFlags::gradients)
645 {
646 const auto *const shape_gradients =
647 shape_data.front().shape_gradients.data();
648 auto *out = fe_eval.begin_gradients();
649 const auto *in = values_dofs_actual;
650
651 for (unsigned int c = 0; c < n_components; c += 3)
652 {
653 if (c + 1 == n_components)
654 apply_matrix_vector_product<evaluate_general,
655 EvaluatorQuantity::value,
656 /*transpose_matrix*/ true,
657 /*add*/ false,
658 /*consider_strides*/ false,
659 Number,
660 Number2,
661 /*n_components*/ 1>(
662 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
663 else if (c + 2 == n_components)
664 apply_matrix_vector_product<evaluate_general,
665 EvaluatorQuantity::value,
666 /*transpose_matrix*/ true,
667 /*add*/ false,
668 /*consider_strides*/ false,
669 Number,
670 Number2,
671 /*n_components*/ 2>(
672 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
673 else
674 apply_matrix_vector_product<evaluate_general,
675 EvaluatorQuantity::value,
676 /*transpose_matrix*/ true,
677 /*add*/ false,
678 /*consider_strides*/ false,
679 Number,
680 Number2,
681 /*n_components*/ 3>(
682 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
683
684 out += 3 * n_q_points * dim;
685 in += 3 * n_dofs;
686 }
687 }
688 }
689
690
691
692 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
693 inline void
694 FEEvaluationImpl<
695 MatrixFreeFunctions::tensor_none,
696 dim,
697 fe_degree,
698 n_q_points_1d,
699 Number>::integrate(const unsigned int n_components,
700 const EvaluationFlags::EvaluationFlags integration_flag,
701 Number *values_dofs_actual,
702 FEEvaluationData<dim, Number, false> &fe_eval,
703 const bool add_into_values_array)
704 {
705 Assert(!(integration_flag & EvaluationFlags::hessians),
706 ExcNotImplemented());
707
708 const std::size_t n_dofs =
709 fe_eval.get_shape_info().dofs_per_component_on_cell;
710 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
711
712 const auto &shape_data = fe_eval.get_shape_info().data;
713
714 using Number2 =
715 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
716
717 if (integration_flag & EvaluationFlags::values)
718 {
719 const auto *const shape_values = shape_data.front().shape_values.data();
720 auto *in = fe_eval.begin_values();
721 auto *out = values_dofs_actual;
722
723 for (unsigned int c = 0; c < n_components; c += 3)
724 {
725 if (add_into_values_array == false)
726 {
727 if (c + 1 == n_components)
728 apply_matrix_vector_product<evaluate_general,
729 EvaluatorQuantity::value,
730 /*transpose_matrix*/ false,
731 /*add*/ false,
732 /*consider_strides*/ false,
733 Number,
734 Number2,
735 /*n_components*/ 1>(
736 shape_values, in, out, n_dofs, n_q_points, 1, 1);
737 else if (c + 2 == n_components)
738 apply_matrix_vector_product<evaluate_general,
739 EvaluatorQuantity::value,
740 /*transpose_matrix*/ false,
741 /*add*/ false,
742 /*consider_strides*/ false,
743 Number,
744 Number2,
745 /*n_components*/ 2>(
746 shape_values, in, out, n_dofs, n_q_points, 1, 1);
747 else
748 apply_matrix_vector_product<evaluate_general,
749 EvaluatorQuantity::value,
750 /*transpose_matrix*/ false,
751 /*add*/ false,
752 /*consider_strides*/ false,
753 Number,
754 Number2,
755 /*n_components*/ 3>(
756 shape_values, in, out, n_dofs, n_q_points, 1, 1);
757 }
758 else
759 {
760 if (c + 1 == n_components)
761 apply_matrix_vector_product<evaluate_general,
762 EvaluatorQuantity::value,
763 /*transpose_matrix*/ false,
764 /*add*/ true,
765 /*consider_strides*/ false,
766 Number,
767 Number2,
768 /*n_components*/ 1>(
769 shape_values, in, out, n_dofs, n_q_points, 1, 1);
770 else if (c + 2 == n_components)
771 apply_matrix_vector_product<evaluate_general,
772 EvaluatorQuantity::value,
773 /*transpose_matrix*/ false,
774 /*add*/ true,
775 /*consider_strides*/ false,
776 Number,
777 Number2,
778 /*n_components*/ 2>(
779 shape_values, in, out, n_dofs, n_q_points, 1, 1);
780 else
781 apply_matrix_vector_product<evaluate_general,
782 EvaluatorQuantity::value,
783 /*transpose_matrix*/ false,
784 /*add*/ true,
785 /*consider_strides*/ false,
786 Number,
787 Number2,
788 /*n_components*/ 3>(
789 shape_values, in, out, n_dofs, n_q_points, 1, 1);
790 }
791 out += 3 * n_dofs;
792 in += 3 * n_q_points;
793 }
794 }
795
796 if (integration_flag & EvaluationFlags::gradients)
797 {
798 const auto *const shape_gradients =
799 shape_data.front().shape_gradients.data();
800 auto *in = fe_eval.begin_gradients();
801 auto *out = values_dofs_actual;
802
803 for (unsigned int c = 0; c < n_components; c += 3)
804 {
805 if (add_into_values_array == false &&
806 !(integration_flag & EvaluationFlags::values))
807 {
808 if (c + 1 == n_components)
809 apply_matrix_vector_product<evaluate_general,
810 EvaluatorQuantity::value,
811 /*transpose_matrix*/ false,
812 /*add*/ false,
813 /*consider_strides*/ false,
814 Number,
815 Number2,
816 /*n_components*/ 1>(
817 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
818 else if (c + 2 == n_components)
819 apply_matrix_vector_product<evaluate_general,
820 EvaluatorQuantity::value,
821 /*transpose_matrix*/ false,
822 /*add*/ false,
823 /*consider_strides*/ false,
824 Number,
825 Number2,
826 /*n_components*/ 2>(
827 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
828 else
829 apply_matrix_vector_product<evaluate_general,
830 EvaluatorQuantity::value,
831 /*transpose_matrix*/ false,
832 /*add*/ false,
833 /*consider_strides*/ false,
834 Number,
835 Number2,
836 /*n_components*/ 3>(
837 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
838 }
839 else
840 {
841 if (c + 1 == n_components)
842 apply_matrix_vector_product<evaluate_general,
843 EvaluatorQuantity::value,
844 /*transpose_matrix*/ false,
845 /*add*/ true,
846 /*consider_strides*/ false,
847 Number,
848 Number2,
849 /*n_components*/ 1>(
850 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
851 else if (c + 2 == n_components)
852 apply_matrix_vector_product<evaluate_general,
853 EvaluatorQuantity::value,
854 /*transpose_matrix*/ false,
855 /*add*/ true,
856 /*consider_strides*/ false,
857 Number,
858 Number2,
859 /*n_components*/ 2>(
860 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
861 else
862 apply_matrix_vector_product<evaluate_general,
863 EvaluatorQuantity::value,
864 /*transpose_matrix*/ false,
865 /*add*/ true,
866 /*consider_strides*/ false,
867 Number,
868 Number2,
869 /*n_components*/ 3>(
870 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
871 }
872 out += 3 * n_dofs;
873 in += 3 * n_q_points * dim;
874 }
875 }
876 }
877
878
879
889 template <EvaluatorVariant variant,
890 EvaluatorQuantity quantity,
891 int dim,
892 int basis_size_1,
893 int basis_size_2>
894 struct FEEvaluationImplBasisChange
895 {
896 static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
897 "The second dimension must not be smaller than the first");
898
921 template <typename Number, typename Number2>
922#ifndef DEBUG
923 DEAL_II_ALWAYS_INLINE
924#endif
925 static void
926 do_forward(const unsigned int n_components,
927 const AlignedVector<Number2> &transformation_matrix,
928 const Number *values_in,
929 Number *values_out,
930 const unsigned int basis_size_1_variable =
931 numbers::invalid_unsigned_int,
932 const unsigned int basis_size_2_variable =
933 numbers::invalid_unsigned_int)
934 {
935 Assert(
936 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
937 ExcMessage("The second dimension must not be smaller than the first"));
938
939 Assert(quantity == EvaluatorQuantity::value, ExcInternalError());
940
941 // we do recursion until dim==1 or dim==2 and we have
942 // basis_size_1==basis_size_2. The latter optimization increases
943 // optimization possibilities for the compiler but does only work for
944 // aliased pointers if the sizes are equal.
945 constexpr int next_dim = (dim == 1 || (dim == 2 && basis_size_1 > 0 &&
946 basis_size_1 == basis_size_2)) ?
947 dim :
948 dim - 1;
949
950 EvaluatorTensorProduct<variant,
951 dim,
952 basis_size_1,
953 (basis_size_1 == 0 ? 0 : basis_size_2),
954 Number,
955 Number2>
956 eval_val(transformation_matrix,
957 {},
958 {},
959 basis_size_1_variable,
960 basis_size_2_variable);
961 const unsigned int np_1 =
962 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
963 const unsigned int np_2 =
964 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
965 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
966 ExcMessage("Cannot transform with 0-point basis"));
967 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
968 ExcMessage("Cannot transform with 0-point basis"));
969
970 // run loop backwards to ensure correctness if values_in aliases with
971 // values_out in case with basis_size_1 < basis_size_2
972 values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
973 values_out =
974 values_out + n_components * Utilities::fixed_power<dim>(np_2);
975 for (unsigned int c = n_components; c != 0; --c)
976 {
977 values_in -= Utilities::fixed_power<dim>(np_1);
978 values_out -= Utilities::fixed_power<dim>(np_2);
979 if (next_dim < dim)
980 for (unsigned int q = np_1; q != 0; --q)
981 FEEvaluationImplBasisChange<variant,
982 quantity,
983 next_dim,
984 basis_size_1,
985 basis_size_2>::
986 do_forward(1,
987 transformation_matrix,
988 values_in +
989 (q - 1) * Utilities::fixed_power<next_dim>(np_1),
990 values_out +
991 (q - 1) * Utilities::fixed_power<next_dim>(np_2),
992 basis_size_1_variable,
993 basis_size_2_variable);
994
995 // the recursion stops if dim==1 or if dim==2 and
996 // basis_size_1==basis_size_2 (the latter is used because the
997 // compiler generates nicer code)
998 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
999 {
1000 eval_val.template values<0, true, false>(values_in, values_out);
1001 eval_val.template values<1, true, false>(values_out, values_out);
1002 }
1003 else if (dim == 1)
1004 eval_val.template values<dim - 1, true, false>(values_in,
1005 values_out);
1006 else
1007 eval_val.template values<dim - 1, true, false>(values_out,
1008 values_out);
1009 }
1010 }
1011
1042 template <typename Number, typename Number2>
1043#ifndef DEBUG
1044 DEAL_II_ALWAYS_INLINE
1045#endif
1046 static void
1047 do_backward(const unsigned int n_components,
1048 const AlignedVector<Number2> &transformation_matrix,
1049 const bool add_into_result,
1050 Number *values_in,
1051 Number *values_out,
1052 const unsigned int basis_size_1_variable =
1053 numbers::invalid_unsigned_int,
1054 const unsigned int basis_size_2_variable =
1055 numbers::invalid_unsigned_int)
1056 {
1057 Assert(
1058 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1059 ExcMessage("The second dimension must not be smaller than the first"));
1060 Assert(add_into_result == false || values_in != values_out,
1061 ExcMessage(
1062 "Input and output cannot alias with each other when "
1063 "adding the result of the basis change to existing data"));
1064
1065 Assert(quantity == EvaluatorQuantity::value ||
1066 quantity == EvaluatorQuantity::hessian,
1067 ExcInternalError());
1068
1069 constexpr int next_dim =
1070 (dim > 2 ||
1071 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1072 dim - 1 :
1073 dim;
1074 EvaluatorTensorProduct<variant,
1075 dim,
1076 basis_size_1,
1077 (basis_size_1 == 0 ? 0 : basis_size_2),
1078 Number,
1079 Number2>
1080 eval_val(transformation_matrix,
1081 transformation_matrix,
1082 transformation_matrix,
1083 basis_size_1_variable,
1084 basis_size_2_variable);
1085 const unsigned int np_1 =
1086 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1087 const unsigned int np_2 =
1088 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1089 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1090 ExcMessage("Cannot transform with 0-point basis"));
1091 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1092 ExcMessage("Cannot transform with 0-point basis"));
1093
1094 for (unsigned int c = 0; c < n_components; ++c)
1095 {
1096 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1097 {
1098 if (quantity == EvaluatorQuantity::value)
1099 eval_val.template values<1, false, false>(values_in, values_in);
1100 else
1101 eval_val.template hessians<1, false, false>(values_in,
1102 values_in);
1103
1104 if (add_into_result)
1105 {
1106 if (quantity == EvaluatorQuantity::value)
1107 eval_val.template values<0, false, true>(values_in,
1108 values_out);
1109 else
1110 eval_val.template hessians<0, false, true>(values_in,
1111 values_out);
1112 }
1113 else
1114 {
1115 if (quantity == EvaluatorQuantity::value)
1116 eval_val.template values<0, false, false>(values_in,
1117 values_out);
1118 else
1119 eval_val.template hessians<0, false, false>(values_in,
1120 values_out);
1121 }
1122 }
1123 else
1124 {
1125 if (dim == 1 && add_into_result)
1126 {
1127 if (quantity == EvaluatorQuantity::value)
1128 eval_val.template values<0, false, true>(values_in,
1129 values_out);
1130 else
1131 eval_val.template hessians<0, false, true>(values_in,
1132 values_out);
1133 }
1134 else if (dim == 1)
1135 {
1136 if (quantity == EvaluatorQuantity::value)
1137 eval_val.template values<0, false, false>(values_in,
1138 values_out);
1139 else
1140 eval_val.template hessians<0, false, false>(values_in,
1141 values_out);
1142 }
1143 else
1144 {
1145 if (quantity == EvaluatorQuantity::value)
1146 eval_val.template values<dim - 1, false, false>(values_in,
1147 values_in);
1148 else
1149 eval_val.template hessians<dim - 1, false, false>(
1150 values_in, values_in);
1151 }
1152 }
1153 if (next_dim < dim)
1154 for (unsigned int q = 0; q < np_1; ++q)
1155 FEEvaluationImplBasisChange<variant,
1156 quantity,
1157 next_dim,
1158 basis_size_1,
1159 basis_size_2>::
1160 do_backward(1,
1161 transformation_matrix,
1162 add_into_result,
1163 values_in +
1164 q * Utilities::fixed_power<next_dim>(np_2),
1165 values_out +
1166 q * Utilities::fixed_power<next_dim>(np_1),
1167 basis_size_1_variable,
1168 basis_size_2_variable);
1169
1170 values_in += Utilities::fixed_power<dim>(np_2);
1171 values_out += Utilities::fixed_power<dim>(np_1);
1172 }
1173 }
1174
1195 template <typename Number, typename Number2>
1196 static void
1197 do_mass(const unsigned int n_components,
1198 const AlignedVector<Number2> &transformation_matrix,
1199 const AlignedVector<Number> &coefficients,
1200 const Number *values_in,
1201 Number *scratch_data,
1202 Number *values_out)
1203 {
1204 constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1205 Number *my_scratch =
1206 basis_size_1 != basis_size_2 ? scratch_data : values_out;
1207
1208 const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1209 Assert(coefficients.size() == size_per_component ||
1210 coefficients.size() == n_components * size_per_component,
1211 ExcDimensionMismatch(coefficients.size(), size_per_component));
1212 const unsigned int stride =
1213 coefficients.size() == size_per_component ? 0 : 1;
1214
1215 for (unsigned int q = basis_size_1; q != 0; --q)
1216 FEEvaluationImplBasisChange<
1217 variant,
1218 EvaluatorQuantity::value,
1219 next_dim,
1220 basis_size_1,
1221 basis_size_2>::do_forward(n_components,
1222 transformation_matrix,
1223 values_in +
1224 (q - 1) *
1225 Utilities::pow(basis_size_1, dim - 1),
1226 my_scratch +
1227 (q - 1) *
1228 Utilities::pow(basis_size_2, dim - 1));
1229 EvaluatorTensorProduct<variant,
1230 dim,
1231 basis_size_1,
1232 basis_size_2,
1233 Number,
1234 Number2>
1235 eval_val(transformation_matrix);
1236 const unsigned int n_inner_blocks =
1237 (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1238 const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1239 for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1240 for (unsigned int c = 0; c < n_components; ++c)
1241 {
1242 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1243 eval_val.template values_one_line<dim - 1, true, false>(
1244 my_scratch + i, my_scratch + i);
1245 for (unsigned int q = 0; q < basis_size_2; ++q)
1246 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1247 my_scratch[i + q * n_blocks + c * size_per_component] *=
1248 coefficients[i + q * n_blocks +
1249 c * stride * size_per_component];
1250 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1251 eval_val.template values_one_line<dim - 1, false, false>(
1252 my_scratch + i, my_scratch + i);
1253 }
1254 for (unsigned int q = 0; q < basis_size_1; ++q)
1255 FEEvaluationImplBasisChange<variant,
1256 EvaluatorQuantity::value,
1257 next_dim,
1258 basis_size_1,
1259 basis_size_2>::
1260 do_backward(n_components,
1261 transformation_matrix,
1262 false,
1263 my_scratch + q * Utilities::pow(basis_size_2, dim - 1),
1264 values_out + q * Utilities::pow(basis_size_1, dim - 1));
1265 }
1266 };
1267
1268
1269
1277 template <int n_points_1d, int dim, typename Number, typename Number2>
1278 inline void
1279 evaluate_gradients_collocation(
1280 const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
1281 const Number *values,
1282 Number *gradients)
1283 {
1284 AssertDimension(shape.shape_gradients_collocation_eo.size(),
1285 (n_points_1d + 1) / 2 * n_points_1d);
1286
1287 EvaluatorTensorProduct<evaluate_evenodd,
1288 dim,
1289 n_points_1d,
1290 n_points_1d,
1291 Number,
1292 Number2>
1293 eval({}, shape.shape_gradients_collocation_eo, {});
1294 EvaluatorTensorProduct<evaluate_evenodd,
1295 2,
1296 n_points_1d,
1297 n_points_1d,
1298 Number,
1299 Number2>
1300 eval_2d({}, shape.shape_gradients_collocation_eo, {});
1301
1302 if (dim == 1)
1303 eval.template gradients<0, true, false>(values, gradients);
1304 else
1305 {
1306 if (dim > 2)
1307 eval.template gradients<2, true, false, dim>(values, gradients + 2);
1308 constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
1309 constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
1310 const Number *in = values + (loop_bound - 1) * n_points_2d;
1311 Number *out = gradients + (loop_bound - 1) * dim * n_points_2d;
1312 for (unsigned int l = 0; l < loop_bound; ++l)
1313 {
1314 eval_2d.template gradients<0, true, false, dim>(in, out);
1315 eval_2d.template gradients<1, true, false, dim>(in, out + 1);
1316 in -= n_points_2d;
1317 out -= dim * n_points_2d;
1318 }
1319 }
1320 }
1321
1322
1323
1331 template <int n_points_1d, int dim, typename Number, typename Number2>
1332 inline void
1333 integrate_gradients_collocation(
1334 const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
1335 Number *values,
1336 const Number *gradients,
1337 const bool add_into_values_array)
1338 {
1339 AssertDimension(shape.shape_gradients_collocation_eo.size(),
1340 (n_points_1d + 1) / 2 * n_points_1d);
1341
1342 EvaluatorTensorProduct<evaluate_evenodd,
1343 dim,
1344 n_points_1d,
1345 n_points_1d,
1346 Number,
1347 Number2>
1348 eval({}, shape.shape_gradients_collocation_eo, {});
1349 EvaluatorTensorProduct<evaluate_evenodd,
1350 2,
1351 n_points_1d,
1352 n_points_1d,
1353 Number,
1354 Number2>
1355 eval_2d({}, shape.shape_gradients_collocation_eo, {});
1356
1357 if (dim == 1)
1358 {
1359 if (add_into_values_array)
1360 eval.template gradients<0, false, true>(gradients, values);
1361 else
1362 eval.template gradients<0, false, false>(gradients, values);
1363 }
1364 else
1365 {
1366 constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
1367 constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
1368
1369 const Number *in = gradients + (loop_bound - 1) * dim * n_points_2d;
1370 Number *out = values + (loop_bound - 1) * n_points_2d;
1371 for (unsigned int l = 0; l < loop_bound; ++l)
1372 {
1373 if (add_into_values_array)
1374 eval_2d.template gradients<0, false, true, dim>(in, out);
1375 else
1376 eval_2d.template gradients<0, false, false, dim>(in, out);
1377 eval_2d.template gradients<1, false, true, dim>(in + 1, out);
1378 in -= dim * n_points_2d;
1379 out -= n_points_2d;
1380 }
1381 }
1382 if (dim > 2)
1383 eval.template gradients<2, false, true, dim>(gradients + 2, values);
1384 }
1385
1386
1387
1394 template <int n_points_1d, int dim, typename Number>
1395 inline void
1396 evaluate_hessians_collocation(const unsigned int n_components,
1397 FEEvaluationData<dim, Number, false> &fe_eval)
1398 {
1399 using Number2 =
1400 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1401
1402 // might have non-symmetric quadrature formula, so use the more
1403 // conservative 'evaluate_general' scheme rather than 'even_odd' as the
1404 // Hessians are not used very often
1405 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data =
1406 fe_eval.get_shape_info().data[0];
1407 AssertDimension(data.shape_gradients_collocation.size(),
1408 data.n_q_points_1d * data.n_q_points_1d);
1409 EvaluatorTensorProduct<evaluate_general,
1410 dim,
1411 n_points_1d,
1412 n_points_1d,
1413 Number,
1414 Number2>
1415 eval({},
1416 data.shape_gradients_collocation.data(),
1417 data.shape_hessians_collocation.data(),
1418 data.n_q_points_1d,
1419 data.n_q_points_1d);
1420
1421 const Number *values = fe_eval.begin_values();
1422 Number *hessians = fe_eval.begin_hessians();
1423 Number *scratch = fe_eval.get_scratch_data().begin();
1424 const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
1425 for (unsigned int comp = 0; comp < n_components; ++comp)
1426 {
1427 // xx derivative
1428 eval.template hessians<0, true, false>(values, hessians);
1429 if (dim > 1)
1430 {
1431 // xy derivative: we might or might not have the gradients already
1432 // computed elsewhere, but we recompute them here since it adds
1433 // only moderate extra work (at most 25%)
1434 eval.template gradients<0, true, false>(values, scratch);
1435 eval.template gradients<1, true, false>(scratch,
1436 hessians + dim * n_points);
1437 // yy derivative
1438 eval.template hessians<1, true, false>(values, hessians + n_points);
1439 }
1440 if (dim > 2)
1441 {
1442 // xz derivative
1443 eval.template gradients<2, true, false>(scratch,
1444 hessians + 4 * n_points);
1445 // yz derivative
1446 eval.template gradients<1, true, false>(values, scratch);
1447 eval.template gradients<2, true, false>(scratch,
1448 hessians + 5 * n_points);
1449 // zz derivative
1450 eval.template hessians<2, true, false>(values,
1451 hessians + 2 * n_points);
1452 }
1453
1454 values += n_points;
1455 hessians += (dim * (dim + 1)) / 2 * n_points;
1456 }
1457 }
1458
1459
1460
1467 template <int n_q_points_1d, int dim, typename Number>
1468 inline void
1469 integrate_hessians_collocation(const unsigned int n_components,
1470 FEEvaluationData<dim, Number, false> &fe_eval,
1471 const bool add_into_values_array)
1472 {
1473 using Number2 =
1474 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1475
1476 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data =
1477 fe_eval.get_shape_info().data[0];
1478 AssertDimension(data.shape_gradients_collocation.size(),
1479 data.n_q_points_1d * data.n_q_points_1d);
1480 EvaluatorTensorProduct<evaluate_general,
1481 dim,
1482 n_q_points_1d,
1483 n_q_points_1d,
1484 Number,
1485 Number2>
1486 eval({},
1487 data.shape_gradients_collocation.data(),
1488 data.shape_hessians_collocation.data(),
1489 data.n_q_points_1d,
1490 data.n_q_points_1d);
1491 Number *values = fe_eval.begin_values();
1492 const Number *hessians = fe_eval.begin_hessians();
1493 Number *scratch = fe_eval.get_scratch_data().begin();
1494 const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
1495
1496 for (unsigned int comp = 0; comp < n_components; ++comp)
1497 {
1498 // xx derivative
1499 if (add_into_values_array == true)
1500 eval.template hessians<0, false, true>(hessians, values);
1501 else
1502 eval.template hessians<0, false, false>(hessians, values);
1503
1504 // yy derivative
1505 if (dim > 1)
1506 eval.template hessians<1, false, true>(hessians + n_points, values);
1507 if (dim > 2)
1508 {
1509 // zz derivative
1510 eval.template hessians<2, false, true>(hessians + 2 * n_points,
1511 values);
1512 // yz derivative
1513 eval.template gradients<2, false, false>(hessians + 5 * n_points,
1514 scratch);
1515 eval.template gradients<1, false, true>(scratch, values);
1516
1517 // xz derivative
1518 eval.template gradients<2, false, false>(hessians + 4 * n_points,
1519 scratch);
1520 }
1521
1522 if (dim > 1)
1523 {
1524 // xy derivative, combined with xz in 3d
1525 eval.template gradients<1, false, (dim > 2)>(hessians +
1526 dim * n_points,
1527 scratch);
1528 eval.template gradients<0, false, true>(scratch, values);
1529 }
1530
1531 values += n_points;
1532 hessians += (dim * (dim + 1)) / 2 * n_points;
1533 }
1534 }
1535
1536
1537
1544 template <int dim, typename Number>
1545 void
1546 evaluate_hessians_slow(const unsigned int n_components,
1547 const Number *values_dofs,
1548 FEEvaluationData<dim, Number, false> &fe_eval)
1549 {
1550 const auto &univariate_shape_data = fe_eval.get_shape_info().data;
1551 using Impl =
1552 FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
1553 using Eval = typename Impl::Eval;
1554 Eval eval0 =
1555 Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
1556 Eval eval1 = Impl::create_evaluator_tensor_product(
1557 &univariate_shape_data[std::min<int>(1,
1558 univariate_shape_data.size() - 1)]);
1559 Eval eval2 = Impl::create_evaluator_tensor_product(
1560 &univariate_shape_data[std::min<int>(2,
1561 univariate_shape_data.size() - 1)]);
1562
1563 const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
1564 Number *tmp1 = fe_eval.get_scratch_data().begin();
1565 Number *tmp2 =
1566 tmp1 + std::max(Utilities::fixed_power<dim>(
1567 univariate_shape_data.front().fe_degree + 1),
1568 Utilities::fixed_power<dim>(
1569 univariate_shape_data.front().n_q_points_1d));
1570 Number *hessians = fe_eval.begin_hessians();
1571
1572 for (unsigned int comp = 0; comp < n_components;
1573 ++comp,
1574 hessians += n_points * dim * (dim + 1) / 2,
1575 values_dofs +=
1576 fe_eval.get_shape_info().dofs_per_component_on_cell)
1577 switch (dim)
1578 {
1579 case 1:
1580 eval0.template hessians<0, true, false>(values_dofs, hessians);
1581 break;
1582 case 2:
1583 // xx derivative
1584 eval0.template hessians<0, true, false>(values_dofs, tmp1);
1585 eval1.template values<1, true, false>(tmp1, hessians);
1586 // xy derivative
1587 eval0.template gradients<0, true, false>(values_dofs, tmp1);
1588 eval1.template gradients<1, true, false>(tmp1,
1589 hessians + 2 * n_points);
1590 // yy derivative
1591 eval0.template values<0, true, false>(values_dofs, tmp1);
1592 eval1.template hessians<1, true, false>(tmp1, hessians + n_points);
1593 break;
1594 case 3:
1595 // xx derivative
1596 eval0.template hessians<0, true, false>(values_dofs, tmp1);
1597 eval1.template values<1, true, false>(tmp1, tmp2);
1598 eval2.template values<2, true, false>(tmp2, hessians);
1599 // xy derivative
1600 eval0.template gradients<0, true, false>(values_dofs, tmp1);
1601 eval1.template gradients<1, true, false>(tmp1, tmp2);
1602 eval2.template values<2, true, false>(tmp2,
1603 hessians + 3 * n_points);
1604 // xz derivative
1605 eval1.template values<1, true, false>(tmp1, tmp2);
1606 eval2.template gradients<2, true, false>(tmp2,
1607 hessians + 4 * n_points);
1608 // yy derivative
1609 eval0.template values<0, true, false>(values_dofs, tmp1);
1610 eval1.template hessians<1, true, false>(tmp1, tmp2);
1611 eval2.template values<2, true, false>(tmp2, hessians + n_points);
1612 // yz derivative
1613 eval1.template gradients<1, true, false>(tmp1, tmp2);
1614 eval2.template gradients<2, true, false>(tmp2,
1615 hessians + 5 * n_points);
1616 // zz derivative
1617 eval1.template values<1, true, false>(tmp1, tmp2);
1618 eval2.template hessians<2, true, false>(tmp2,
1619 hessians + 2 * n_points);
1620 break;
1621
1622 default:
1623 Assert(false,
1624 ExcNotImplemented(
1625 "Only 1d, 2d and 3d implemented for Hessian"));
1626 }
1627 }
1628
1629
1630
1638 template <int dim, typename Number>
1639 void
1640 integrate_hessians_slow(const unsigned int n_components,
1641 const FEEvaluationData<dim, Number, false> &fe_eval,
1642 Number *values_dofs,
1643 const bool add_into_values_array)
1644 {
1645 const auto &univariate_shape_data = fe_eval.get_shape_info().data;
1646 using Impl =
1647 FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
1648 using Eval = typename Impl::Eval;
1649 Eval eval0 =
1650 Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
1651 Eval eval1 = Impl::create_evaluator_tensor_product(
1652 &univariate_shape_data[std::min<int>(1,
1653 univariate_shape_data.size() - 1)]);
1654 Eval eval2 = Impl::create_evaluator_tensor_product(
1655 &univariate_shape_data[std::min<int>(2,
1656 univariate_shape_data.size() - 1)]);
1657
1658 const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
1659 Number *tmp1 = fe_eval.get_scratch_data().begin();
1660 Number *tmp2 =
1661 tmp1 + std::max(Utilities::fixed_power<dim>(
1662 univariate_shape_data.front().fe_degree + 1),
1663 Utilities::fixed_power<dim>(
1664 univariate_shape_data.front().n_q_points_1d));
1665 const Number *hessians = fe_eval.begin_hessians();
1666
1667 for (unsigned int comp = 0; comp < n_components;
1668 ++comp,
1669 hessians += n_points * dim * (dim + 1) / 2,
1670 values_dofs +=
1671 fe_eval.get_shape_info().dofs_per_component_on_cell)
1672 switch (dim)
1673 {
1674 case 1:
1675 if (add_into_values_array)
1676 eval0.template hessians<0, false, true>(hessians, values_dofs);
1677 else
1678 eval0.template hessians<0, false, false>(hessians, values_dofs);
1679 break;
1680 case 2:
1681 // xx derivative
1682 eval1.template values<1, false, false>(hessians, tmp1);
1683 if (add_into_values_array)
1684 eval0.template hessians<0, false, true>(tmp1, values_dofs);
1685 else
1686 eval0.template hessians<0, false, false>(tmp1, values_dofs);
1687
1688 // xy derivative
1689 eval1.template gradients<1, false, false>(hessians + 2 * n_points,
1690 tmp1);
1691 eval0.template gradients<0, false, true>(tmp1, values_dofs);
1692 // yy derivative
1693 eval1.template hessians<1, false, false>(hessians + n_points, tmp1);
1694 eval0.template values<0, false, true>(tmp1, values_dofs);
1695 break;
1696 case 3:
1697 // xx derivative
1698 eval2.template values<2, false, false>(hessians, tmp1);
1699 eval1.template values<1, false, false>(tmp1, tmp2);
1700
1701 if (add_into_values_array)
1702 eval0.template hessians<0, false, true>(tmp2, values_dofs);
1703 else
1704 eval0.template hessians<0, false, false>(tmp2, values_dofs);
1705
1706 // xy derivative
1707 eval2.template values<2, false, false>(hessians + 3 * n_points,
1708 tmp1);
1709 eval1.template gradients<1, false, false>(tmp1, tmp2);
1710 // xz derivative
1711 eval2.template gradients<2, false, false>(hessians + 4 * n_points,
1712 tmp1);
1713 eval1.template values<1, false, true>(tmp1, tmp2);
1714 eval1.template values<0, false, true>(tmp2, values_dofs);
1715
1716 // yy derivative
1717 eval2.template values<2, false, false>(hessians + n_points, tmp1);
1718 eval1.template hessians<1, false, false>(tmp1, tmp2);
1719
1720 // yz derivative
1721 eval2.template gradients<2, false, false>(hessians + 5 * n_points,
1722 tmp1);
1723 eval1.template gradients<1, false, true>(tmp1, tmp2);
1724
1725 // zz derivative
1726 eval2.template hessians<2, false, false>(hessians + 2 * n_points,
1727 tmp1);
1728 eval1.template values<1, false, true>(tmp1, tmp2);
1729 eval0.template values<0, false, true>(tmp2, values_dofs);
1730 break;
1731
1732 default:
1733 Assert(false,
1734 ExcNotImplemented(
1735 "Only 1d, 2d and 3d implemented for Hessian"));
1736 }
1737 }
1738
1739
1740
1753 template <int dim, int fe_degree, typename Number>
1754 struct FEEvaluationImplCollocation
1755 {
1756 using Number2 =
1757 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1758 using Eval = EvaluatorTensorProduct<evaluate_evenodd,
1759 dim,
1760 fe_degree + 1,
1761 fe_degree + 1,
1762 Number,
1763 Number2>;
1764
1765 static void
1766 evaluate(const unsigned int n_components,
1767 const EvaluationFlags::EvaluationFlags evaluation_flag,
1768 const Number *values_dofs,
1769 FEEvaluationData<dim, Number, false> &fe_eval)
1770 {
1771 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1772
1773 for (unsigned int c = 0; c < n_components; ++c)
1774 {
1775 if ((evaluation_flag & EvaluationFlags::values) != 0u)
1776 for (unsigned int i = 0; i < n_points; ++i)
1777 fe_eval.begin_values()[n_points * c + i] =
1778 values_dofs[n_points * c + i];
1779
1780 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1781 evaluate_gradients_collocation<fe_degree + 1, dim>(
1782 fe_eval.get_shape_info().data.front(),
1783 values_dofs + c * n_points,
1784 fe_eval.begin_gradients() + c * dim * n_points);
1785 }
1786 }
1787
1788 static void
1789 integrate(const unsigned int n_components,
1790 const EvaluationFlags::EvaluationFlags integration_flag,
1791 Number *values_dofs,
1792 FEEvaluationData<dim, Number, false> &fe_eval,
1793 const bool add_into_values_array)
1794 {
1795 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1796
1797 for (unsigned int c = 0; c < n_components; ++c)
1798 {
1799 if ((integration_flag & EvaluationFlags::values) != 0u)
1800 {
1801 if (add_into_values_array)
1802 for (unsigned int i = 0; i < n_points; ++i)
1803 values_dofs[n_points * c + i] +=
1804 fe_eval.begin_values()[n_points * c + i];
1805 else
1806 for (unsigned int i = 0; i < n_points; ++i)
1807 values_dofs[n_points * c + i] =
1808 fe_eval.begin_values()[n_points * c + i];
1809 }
1810
1811 if ((integration_flag & EvaluationFlags::gradients) != 0u)
1812 integrate_gradients_collocation<fe_degree + 1, dim>(
1813 fe_eval.get_shape_info().data.front(),
1814 values_dofs + c * n_points,
1815 fe_eval.begin_gradients() + c * dim * n_points,
1816 add_into_values_array ||
1817 ((integration_flag & EvaluationFlags::values) != 0u));
1818 }
1819 }
1820 };
1821
1822
1823
1834 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1835 struct FEEvaluationImplTransformToCollocation
1836 {
1837 static void
1838 evaluate(const unsigned int n_components,
1839 const EvaluationFlags::EvaluationFlags evaluation_flag,
1840 const Number *values_dofs,
1841 FEEvaluationData<dim, Number, false> &fe_eval)
1842 {
1843 const auto &shape_data = fe_eval.get_shape_info().data.front();
1844
1845 Assert(n_q_points_1d > fe_degree,
1846 ExcMessage("You lose information when going to a collocation "
1847 "space of lower degree, so the evaluation results "
1848 "would be wrong. Thus, this class does not permit "
1849 "the chosen operation."));
1850 constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
1851 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
1852
1853 for (unsigned int c = 0; c < n_components; ++c)
1854 {
1855 FEEvaluationImplBasisChange<
1856 evaluate_evenodd,
1857 EvaluatorQuantity::value,
1858 dim,
1859 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1860 n_q_points_1d>::do_forward(1,
1861 shape_data.shape_values_eo,
1862 values_dofs + c * n_dofs,
1863 fe_eval.begin_values() + c * n_q_points);
1864
1865 // apply derivatives in the collocation space
1866 if (evaluation_flag & EvaluationFlags::gradients)
1867 evaluate_gradients_collocation<n_q_points_1d, dim>(
1868 shape_data,
1869 fe_eval.begin_values() + c * n_q_points,
1870 fe_eval.begin_gradients() + c * dim * n_q_points);
1871 }
1872 }
1873
1874 static void
1875 integrate(const unsigned int n_components,
1876 const EvaluationFlags::EvaluationFlags integration_flag,
1877 Number *values_dofs,
1878 FEEvaluationData<dim, Number, false> &fe_eval,
1879 const bool add_into_values_array)
1880 {
1881 const auto &shape_data = fe_eval.get_shape_info().data.front();
1882
1883 Assert(n_q_points_1d > fe_degree,
1884 ExcMessage("You lose information when going to a collocation "
1885 "space of lower degree, so the evaluation results "
1886 "would be wrong. Thus, this class does not permit "
1887 "the chosen operation."));
1888 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
1889
1890 for (unsigned int c = 0; c < n_components; ++c)
1891 {
1892 // apply derivatives in collocation space
1893 if (integration_flag & EvaluationFlags::gradients)
1894 integrate_gradients_collocation<n_q_points_1d, dim>(
1895 shape_data,
1896 fe_eval.begin_values() + c * n_q_points,
1897 fe_eval.begin_gradients() + c * dim * n_q_points,
1898 /*add_into_values_array=*/
1899 integration_flag & EvaluationFlags::values);
1900
1901 // transform back to the original space
1902 FEEvaluationImplBasisChange<
1903 evaluate_evenodd,
1904 EvaluatorQuantity::value,
1905 dim,
1906 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1907 n_q_points_1d>::do_backward(1,
1908 shape_data.shape_values_eo,
1909 add_into_values_array,
1910 fe_eval.begin_values() + c * n_q_points,
1911 values_dofs +
1912 c *
1913 Utilities::pow(fe_degree + 1, dim));
1914 }
1915 }
1916 };
1917
1918
1919
1924 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1925 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
1926 dim,
1927 fe_degree,
1928 n_q_points_1d,
1929 Number>
1930 {
1931 using Number2 =
1932 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1933
1934 template <bool integrate>
1935 static void
1936 evaluate_or_integrate(
1937 const EvaluationFlags::EvaluationFlags evaluation_flag,
1938 Number *values_dofs_actual,
1939 FEEvaluationData<dim, Number, false> &fe_eval,
1940 const bool add_into_values_array = false);
1941 };
1942
1943
1944
1945 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1946 template <bool integrate>
1947 inline void
1948 FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
1949 dim,
1950 fe_degree,
1951 n_q_points_1d,
1952 Number>::
1953 evaluate_or_integrate(
1954 const EvaluationFlags::EvaluationFlags evaluation_flag,
1955 Number *values_dofs,
1956 FEEvaluationData<dim, Number, false> &fe_eval,
1957 const bool add)
1958 {
1959 Assert(dim == 2 || dim == 3,
1960 ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas "
1961 "evaluation/integration"));
1962
1963 if (evaluation_flag == EvaluationFlags::nothing)
1964 return;
1965
1966 AssertDimension(fe_eval.get_shape_info().data.size(), 2);
1967 AssertDimension(n_q_points_1d,
1968 fe_eval.get_shape_info().data[0].n_q_points_1d);
1969 AssertDimension(n_q_points_1d,
1970 fe_eval.get_shape_info().data[1].n_q_points_1d);
1971 AssertDimension(fe_degree, fe_eval.get_shape_info().data[0].fe_degree);
1972 AssertDimension(fe_degree, fe_eval.get_shape_info().data[1].fe_degree + 1);
1973
1974 const auto &shape_data = fe_eval.get_shape_info().data;
1975 const unsigned int dofs_per_component =
1976 Utilities::pow(fe_degree, dim - 1) * (fe_degree + 1);
1977 const unsigned int n_points = Utilities::pow(n_q_points_1d, dim);
1978 Number *gradients = fe_eval.begin_gradients();
1979 Number *values = fe_eval.begin_values();
1980
1981 if (integrate)
1982 {
1983 EvaluatorTensorProductAnisotropic<dim, fe_degree, n_q_points_1d, false>
1984 eval;
1985
1986 const bool do_values = evaluation_flag & EvaluationFlags::values;
1987 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1988 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1989 values,
1990 gradients,
1991 do_values);
1992 if constexpr (dim > 2)
1993 eval.template tangential<2, 0>(shape_data[1], values, values);
1994 eval.template tangential<1, 0>(shape_data[1], values, values);
1995 eval.template normal<0>(shape_data[0], values, values_dofs, add);
1996
1997 values += n_points;
1998 gradients += n_points * dim;
1999 values_dofs += dofs_per_component;
2000
2001 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2002 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2003 values,
2004 gradients,
2005 do_values);
2006 if constexpr (dim > 2)
2007 eval.template tangential<2, 1>(shape_data[1], values, values);
2008 eval.template tangential<0, 1>(shape_data[1], values, values);
2009 eval.template normal<1>(shape_data[0], values, values_dofs, add);
2010
2011 if constexpr (dim > 2)
2012 {
2013 values += n_points;
2014 gradients += n_points * dim;
2015 values_dofs += dofs_per_component;
2016
2017 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2018 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2019 values,
2020 gradients,
2021 do_values);
2022 eval.template tangential<1, 2>(shape_data[1], values, values);
2023 eval.template tangential<0, 2>(shape_data[1], values, values);
2024 eval.template normal<2>(shape_data[0], values, values_dofs, add);
2025 }
2026 }
2027 else
2028 {
2029 EvaluatorTensorProductAnisotropic<dim, fe_degree, n_q_points_1d, true>
2030 eval;
2031 eval.template normal<0>(shape_data[0], values_dofs, values);
2032 eval.template tangential<1, 0>(shape_data[1], values, values);
2033 if constexpr (dim > 2)
2034 eval.template tangential<2, 0>(shape_data[1], values, values);
2035 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2036 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2037 values,
2038 gradients);
2039
2040 values += n_points;
2041 gradients += n_points * dim;
2042 values_dofs += dofs_per_component;
2043
2044 eval.template normal<1>(shape_data[0], values_dofs, values);
2045 eval.template tangential<0, 1>(shape_data[1], values, values);
2046 if constexpr (dim > 2)
2047 eval.template tangential<2, 1>(shape_data[1], values, values);
2048 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2049 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2050 values,
2051 gradients);
2052
2053 if constexpr (dim > 2)
2054 {
2055 values += n_points;
2056 gradients += n_points * dim;
2057 values_dofs += dofs_per_component;
2058
2059 eval.template normal<2>(shape_data[0], values_dofs, values);
2060 eval.template tangential<0, 2>(shape_data[1], values, values);
2061 eval.template tangential<1, 2>(shape_data[1], values, values);
2062 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2063 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2064 values,
2065 gradients);
2066 }
2067 }
2068 }
2069
2070
2071
2076 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2077 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_nedelec,
2078 dim,
2079 fe_degree,
2080 n_q_points_1d,
2081 Number>
2082 {
2083 using Number2 =
2084 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2085
2086 template <bool integrate>
2087 static void
2088 evaluate_or_integrate(
2089 const EvaluationFlags::EvaluationFlags evaluation_flag,
2090 Number *values_dofs_actual,
2091 FEEvaluationData<dim, Number, false> &fe_eval,
2092 const bool add_into_values_array = false);
2093 };
2094
2095
2096
2097 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2098 template <bool integrate>
2099 inline void
2100 FEEvaluationImpl<MatrixFreeFunctions::tensor_nedelec,
2101 dim,
2102 fe_degree,
2103 n_q_points_1d,
2104 Number>::
2105 evaluate_or_integrate(
2106 const EvaluationFlags::EvaluationFlags evaluation_flag,
2107 Number *values_dofs,
2108 FEEvaluationData<dim, Number, false> &fe_eval,
2109 const bool add)
2110 {
2111 Assert(dim == 2 || dim == 3,
2112 ExcMessage("Only dim = 2,3 implemented for Nedelec "
2113 "evaluation/integration"));
2114
2115 if (evaluation_flag == EvaluationFlags::nothing)
2116 return;
2117
2118 AssertDimension(fe_eval.get_shape_info().data.size(), 2);
2119 AssertDimension(n_q_points_1d,
2120 fe_eval.get_shape_info().data[0].n_q_points_1d);
2121 AssertDimension(n_q_points_1d,
2122 fe_eval.get_shape_info().data[1].n_q_points_1d);
2123 AssertDimension(fe_degree, fe_eval.get_shape_info().data[0].fe_degree);
2124 AssertDimension(fe_degree, fe_eval.get_shape_info().data[1].fe_degree - 1);
2125
2126
2127 const auto &shape_data = fe_eval.get_shape_info().data;
2128 const unsigned int dofs_per_component =
2129 Utilities::pow(fe_degree + 2, dim - 1) * (fe_degree + 1);
2130 const unsigned int n_points = Utilities::pow(n_q_points_1d, dim);
2131 Number *gradients = fe_eval.begin_gradients();
2132 Number *values = fe_eval.begin_values();
2133
2134 if (integrate)
2135 {
2136 EvaluatorTensorProductAnisotropic<dim,
2137 fe_degree,
2138 n_q_points_1d,
2139 false,
2140 MatrixFreeFunctions::tensor_nedelec>
2141 eval;
2142
2143 const bool do_values = evaluation_flag & EvaluationFlags::values;
2144 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2145 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2146 values,
2147 gradients,
2148 do_values);
2149
2150 if constexpr (dim > 2)
2151 eval.template tangential<2, 0>(shape_data[1], values, values);
2152 eval.template tangential<1, 0>(shape_data[1], values, values);
2153 eval.template normal<0>(shape_data[0], values, values_dofs, add);
2154
2155
2156
2157 values += n_points;
2158 gradients += n_points * dim;
2159 values_dofs += dofs_per_component;
2160
2161 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2162 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2163 values,
2164 gradients,
2165 do_values);
2166
2167 if constexpr (dim > 2)
2168 eval.template tangential<2, 1>(shape_data[1], values, values);
2169 eval.template tangential<0, 1>(shape_data[1], values, values);
2170 eval.template normal<1>(shape_data[0], values, values_dofs, add);
2171
2172 if constexpr (dim > 2)
2173 {
2174 values += n_points;
2175 gradients += n_points * dim;
2176 values_dofs += dofs_per_component;
2177
2178 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2179 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2180 values,
2181 gradients,
2182 do_values);
2183
2184 eval.template tangential<1, 2>(shape_data[1], values, values);
2185 eval.template tangential<0, 2>(shape_data[1], values, values);
2186 eval.template normal<2>(shape_data[0], values, values_dofs, add);
2187 }
2188 }
2189 else
2190 {
2191 EvaluatorTensorProductAnisotropic<dim,
2192 fe_degree,
2193 n_q_points_1d,
2194 true,
2195 MatrixFreeFunctions::tensor_nedelec>
2196 eval;
2197 eval.template normal<0>(shape_data[0], values_dofs, values);
2198 eval.template tangential<1, 0>(shape_data[1], values, values);
2199 if constexpr (dim > 2)
2200 eval.template tangential<2, 0>(shape_data[1], values, values);
2201 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2202 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2203 values,
2204 gradients);
2205
2206 values += n_points;
2207 gradients += n_points * dim;
2208 values_dofs += dofs_per_component;
2209
2210 eval.template normal<1>(shape_data[0], values_dofs, values);
2211 eval.template tangential<0, 1>(shape_data[1], values, values);
2212 if constexpr (dim > 2)
2213 eval.template tangential<2, 1>(shape_data[1], values, values);
2214 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2215 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2216 values,
2217 gradients);
2218
2219 if constexpr (dim > 2)
2220 {
2221 values += n_points;
2222 gradients += n_points * dim;
2223 values_dofs += dofs_per_component;
2224
2225 eval.template normal<2>(shape_data[0], values_dofs, values);
2226 eval.template tangential<0, 2>(shape_data[1], values, values);
2227 eval.template tangential<1, 2>(shape_data[1], values, values);
2228 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2229 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
2230 values,
2231 gradients);
2232 }
2233 }
2234 }
2235
2236
2237
2253 template <int dim, typename Number, bool do_integrate>
2254 struct FEEvaluationImplSelector
2255 {
2256 template <int fe_degree, int n_q_points_1d, typename OtherNumber>
2257 static bool
2258 run(const unsigned int n_components,
2259 const EvaluationFlags::EvaluationFlags evaluation_flag,
2260 OtherNumber *values_dofs,
2261 FEEvaluationData<dim, Number, false> &fe_eval,
2262 const bool sum_into_values_array_in = false)
2263 {
2264 // `OtherNumber` is either `const Number` (evaluate()) or `Number`
2265 // (integrate())
2266 static_assert(std::is_same_v<Number, std::remove_const_t<OtherNumber>>,
2267 "Type of Number and of OtherNumber do not match.");
2268
2269 const auto element_type = fe_eval.get_shape_info().element_type;
2270 using ElementType = MatrixFreeFunctions::ElementType;
2271
2272 Assert(fe_eval.get_shape_info().data.size() == 1 ||
2273 (fe_eval.get_shape_info().data.size() == dim &&
2274 element_type == ElementType::tensor_general) ||
2275 element_type == ElementType::tensor_raviart_thomas ||
2276 element_type == ElementType::tensor_nedelec,
2277 ExcNotImplemented());
2278
2279 EvaluationFlags::EvaluationFlags actual_flag = evaluation_flag;
2280 bool sum_into_values_array = sum_into_values_array_in;
2281 if (evaluation_flag & EvaluationFlags::hessians)
2282 {
2283 actual_flag |= EvaluationFlags::values;
2284 Assert(element_type != MatrixFreeFunctions::tensor_none,
2285 ExcNotImplemented());
2286 if constexpr (do_integrate)
2287 {
2288 if (fe_eval.get_shape_info().data[0].fe_degree <
2289 fe_eval.get_shape_info().data[0].n_q_points_1d)
2290 integrate_hessians_collocation<n_q_points_1d>(
2291 n_components,
2292 fe_eval,
2293 evaluation_flag & EvaluationFlags::values);
2294 else
2295 {
2296 integrate_hessians_slow(n_components,
2297 fe_eval,
2298 values_dofs,
2299 sum_into_values_array);
2300 sum_into_values_array = true;
2301 }
2302 }
2303 }
2304
2305 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2306 element_type == ElementType::tensor_symmetric_collocation)
2307 {
2308 evaluate_or_integrate<
2309 FEEvaluationImplCollocation<dim, fe_degree, Number>>(
2310 n_components,
2311 actual_flag,
2312 values_dofs,
2313 fe_eval,
2314 sum_into_values_array);
2315 }
2316 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2317 // shape_info.h for more details
2318 else if (fe_degree >= 0 &&
2319 use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2320 element_type <= ElementType::tensor_symmetric)
2321 {
2322 evaluate_or_integrate<
2323 FEEvaluationImplTransformToCollocation<dim,
2324 fe_degree,
2325 n_q_points_1d,
2326 Number>>(
2327 n_components,
2328 actual_flag,
2329 values_dofs,
2330 fe_eval,
2331 sum_into_values_array);
2332 }
2333 else if (fe_degree >= 0 &&
2334 element_type <= ElementType::tensor_symmetric_no_collocation)
2335 {
2336 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_symmetric,
2337 dim,
2338 fe_degree,
2339 n_q_points_1d,
2340 Number>>(
2341 n_components,
2342 actual_flag,
2343 values_dofs,
2344 fe_eval,
2345 sum_into_values_array);
2346 }
2347 else if (element_type == ElementType::tensor_none)
2348 {
2349 evaluate_or_integrate<
2350 FEEvaluationImpl<ElementType::tensor_none, dim, -1, 0, Number>>(
2351 n_components,
2352 actual_flag,
2353 values_dofs,
2354 fe_eval,
2355 sum_into_values_array);
2356 }
2357 else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2358 {
2359 evaluate_or_integrate<
2360 FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
2361 dim,
2362 fe_degree,
2363 n_q_points_1d,
2364 Number>>(n_components,
2365 actual_flag,
2366 values_dofs,
2367 fe_eval,
2368 sum_into_values_array);
2369 }
2370 else if (element_type == ElementType::truncated_tensor)
2371 {
2372 evaluate_or_integrate<FEEvaluationImpl<ElementType::truncated_tensor,
2373 dim,
2374 fe_degree,
2375 n_q_points_1d,
2376 Number>>(
2377 n_components,
2378 actual_flag,
2379 values_dofs,
2380 fe_eval,
2381 sum_into_values_array);
2382 }
2383 else if (element_type == ElementType::tensor_raviart_thomas)
2384 {
2385 if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1)
2386 {
2387 FEEvaluationImpl<ElementType::tensor_raviart_thomas,
2388 dim,
2389 fe_degree,
2390 n_q_points_1d,
2391 Number>::
2392 template evaluate_or_integrate<do_integrate>(
2393 actual_flag,
2394 const_cast<Number *>(values_dofs),
2395 fe_eval,
2396 sum_into_values_array);
2397 }
2398 else
2399 {
2400 AssertThrow(false,
2401 ExcNotImplemented(
2402 "Raviart-Thomas currently only possible "
2403 "in 2d/3d and with templated degree for "
2404 "requested fe_degree and n_q_points_1d. "
2405 "Ensure that the highest used degree for "
2406 "Raviart-Thomas, fe_degree+1=" +
2407 std::to_string(
2408 fe_eval.get_shape_info().data[0].fe_degree) +
2409 ", does not exceed FE_EVAL_FACTORY_DEGREE_MAX."));
2410 }
2411 }
2412 else if (element_type == ElementType::tensor_nedelec)
2413 {
2414 if constexpr (fe_degree >= 0 && n_q_points_1d > 0 && dim > 1)
2415 {
2416 FEEvaluationImpl<ElementType::tensor_nedelec,
2417 dim,
2418 fe_degree,
2419 n_q_points_1d,
2420 Number>::
2421 template evaluate_or_integrate<do_integrate>(
2422 actual_flag,
2423 const_cast<Number *>(values_dofs),
2424 fe_eval,
2425 sum_into_values_array);
2426 }
2427 else
2428 {
2429 Assert(false,
2430 ExcNotImplemented("Nedelec currently only possible "
2431 "in 2d/3d and with templated degree"));
2432 }
2433 }
2434 else
2435 {
2436 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_general,
2437 dim,
2438 fe_degree,
2439 n_q_points_1d,
2440 Number>>(
2441 n_components,
2442 actual_flag,
2443 values_dofs,
2444 fe_eval,
2445 sum_into_values_array);
2446 }
2447
2448 if ((evaluation_flag & EvaluationFlags::hessians) && !do_integrate)
2449 {
2450 Assert(element_type != MatrixFreeFunctions::tensor_none,
2451 ExcNotImplemented());
2452 if (fe_eval.get_shape_info().data[0].fe_degree <
2453 fe_eval.get_shape_info().data[0].n_q_points_1d)
2454 evaluate_hessians_collocation<n_q_points_1d>(n_components, fe_eval);
2455 else
2456 evaluate_hessians_slow(n_components, values_dofs, fe_eval);
2457 }
2458
2459 return false;
2460 }
2461
2462 private:
2463 template <typename T>
2464 static void
2465 evaluate_or_integrate(
2466 const unsigned int n_components,
2467 const EvaluationFlags::EvaluationFlags evaluation_flag,
2468 const Number *values_dofs,
2469 FEEvaluationData<dim, Number, false> &fe_eval,
2470 const bool sum_into_values_array,
2471 std::bool_constant<false>)
2472 {
2473 (void)sum_into_values_array;
2474
2475 T::evaluate(n_components, evaluation_flag, values_dofs, fe_eval);
2476 }
2477
2478 template <typename T>
2479 static void
2480 evaluate_or_integrate(
2481 const unsigned int n_components,
2482 const EvaluationFlags::EvaluationFlags evaluation_flag,
2483 Number *values_dofs,
2484 FEEvaluationData<dim, Number, false> &fe_eval,
2485 const bool sum_into_values_array,
2486 std::bool_constant<true>)
2487 {
2488 T::integrate(n_components,
2489 evaluation_flag,
2490 values_dofs,
2491 fe_eval,
2492 sum_into_values_array);
2493 }
2494
2495 template <typename T, typename OtherNumber>
2496 static void
2497 evaluate_or_integrate(
2498 const unsigned int n_components,
2499 const EvaluationFlags::EvaluationFlags evaluation_flag,
2500 OtherNumber *values_dofs,
2501 FEEvaluationData<dim, Number, false> &fe_eval,
2502 const bool sum_into_values_array)
2503 {
2504 evaluate_or_integrate<T>(n_components,
2505 evaluation_flag,
2506 values_dofs,
2507 fe_eval,
2508 sum_into_values_array,
2509 std::bool_constant<do_integrate>());
2510 }
2511 };
2512
2513
2514
2519 template <int dim, typename Number>
2520 struct CellwiseInverseMassMatrixImplBasic
2521 {
2522 using Number2 =
2523 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2524
2525 template <int fe_degree, int = 0>
2526 static bool
2527 run(const unsigned int n_components,
2528 const FEEvaluationData<dim, Number, false> &fe_eval,
2529 const Number *in_array,
2530 Number *out_array)
2531 {
2532 const unsigned int given_degree =
2533 (fe_degree > -1) ? fe_degree :
2534 fe_eval.get_shape_info().data.front().fe_degree;
2535
2536 const unsigned int dofs_per_component =
2537 Utilities::pow(given_degree + 1, dim);
2538
2539 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
2540 Assert(fe_eval.get_shape_info().element_type <=
2541 MatrixFreeFunctions::tensor_symmetric_no_collocation,
2542 ExcNotImplemented());
2543
2544 EvaluatorTensorProduct<evaluate_evenodd,
2545 dim,
2546 fe_degree + 1,
2547 fe_degree + 1,
2548 Number,
2549 Number2>
2550 evaluator({},
2551 {},
2552 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
2553 given_degree + 1,
2554 given_degree + 1);
2555
2556 for (unsigned int d = 0; d < n_components; ++d)
2557 {
2558 const Number *in = in_array + d * dofs_per_component;
2559 Number *out = out_array + d * dofs_per_component;
2560 // Need to select 'apply' method with hessian slot because values
2561 // assume symmetries that do not exist in the inverse shapes
2562 evaluator.template hessians<0, true, false>(in, out);
2563 if (dim > 1)
2564 evaluator.template hessians<1, true, false>(out, out);
2565 if (dim > 2)
2566 evaluator.template hessians<2, true, false>(out, out);
2567 }
2568 for (unsigned int q = 0; q < dofs_per_component; ++q)
2569 {
2570 const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
2571 for (unsigned int d = 0; d < n_components; ++d)
2572 out_array[q + d * dofs_per_component] *= inverse_JxW_q;
2573 }
2574 for (unsigned int d = 0; d < n_components; ++d)
2575 {
2576 Number *out = out_array + d * dofs_per_component;
2577 if (dim > 2)
2578 evaluator.template hessians<2, false, false>(out, out);
2579 if (dim > 1)
2580 evaluator.template hessians<1, false, false>(out, out);
2581 evaluator.template hessians<0, false, false>(out, out);
2582 }
2583 return false;
2584 }
2585 };
2586
2587
2588
2595 template <int dim, typename Number>
2596 struct CellwiseInverseMassMatrixImplFlexible
2597 {
2598 using Number2 =
2599 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2600
2601 template <int fe_degree, int = 0>
2602 static bool
2603 run(const unsigned int n_desired_components,
2604 const FEEvaluationData<dim, Number, false> &fe_eval,
2605 const ArrayView<const Number> &inverse_coefficients,
2606 const bool dyadic_coefficients,
2607 const Number *in_array,
2608 Number *out_array)
2609 {
2610 const unsigned int given_degree =
2611 (fe_degree > -1) ? fe_degree :
2612 fe_eval.get_shape_info().data.front().fe_degree;
2613
2614 const unsigned int dofs_per_component =
2615 Utilities::pow(given_degree + 1, dim);
2616
2617 Assert(inverse_coefficients.size() > 0 &&
2618 inverse_coefficients.size() % dofs_per_component == 0,
2619 ExcMessage(
2620 "Expected diagonal to be a multiple of scalar dof per cells"));
2621
2622 if (!dyadic_coefficients)
2623 {
2624 if (inverse_coefficients.size() != dofs_per_component)
2625 AssertDimension(n_desired_components * dofs_per_component,
2626 inverse_coefficients.size());
2627 }
2628 else
2629 {
2630 AssertDimension(n_desired_components * n_desired_components *
2631 dofs_per_component,
2632 inverse_coefficients.size());
2633 }
2634
2635 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
2636 Assert(fe_eval.get_shape_info().element_type <=
2637 MatrixFreeFunctions::tensor_symmetric_no_collocation,
2638 ExcNotImplemented());
2639
2640 EvaluatorTensorProduct<evaluate_evenodd,
2641 dim,
2642 fe_degree + 1,
2643 fe_degree + 1,
2644 Number,
2645 Number2>
2646 evaluator({},
2647 {},
2648 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
2649 given_degree + 1,
2650 given_degree + 1);
2651
2652 const Number *in = in_array;
2653 Number *out = out_array;
2654
2655 const Number *inv_coefficient = inverse_coefficients.data();
2656
2657 const unsigned int shift_coefficient =
2658 inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
2659 0;
2660
2661 const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
2662 const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
2663
2664 for (unsigned int d = 0; d < n_comp_outer; ++d)
2665 {
2666 for (unsigned int di = 0; di < n_comp_inner; ++di)
2667 {
2668 const Number *in_ = in + di * dofs_per_component;
2669 Number *out_ = out + di * dofs_per_component;
2670 evaluator.template hessians<0, true, false>(in_, out_);
2671 if (dim > 1)
2672 evaluator.template hessians<1, true, false>(out_, out_);
2673 if (dim > 2)
2674 evaluator.template hessians<2, true, false>(out_, out_);
2675 }
2676 if (dyadic_coefficients)
2677 {
2678 const auto n_coeff_components =
2679 n_desired_components * n_desired_components;
2680 if (n_desired_components == dim)
2681 {
2682 for (unsigned int q = 0; q < dofs_per_component; ++q)
2683 vmult<dim>(&inv_coefficient[q * n_coeff_components],
2684 &in[q],
2685 &out[q],
2686 dofs_per_component);
2687 }
2688 else
2689 {
2690 for (unsigned int q = 0; q < dofs_per_component; ++q)
2691 vmult<-1>(&inv_coefficient[q * n_coeff_components],
2692 &in[q],
2693 &out[q],
2694 dofs_per_component,
2695 n_desired_components);
2696 }
2697 }
2698 else
2699 for (unsigned int q = 0; q < dofs_per_component; ++q)
2700 out[q] *= inv_coefficient[q];
2701
2702 for (unsigned int di = 0; di < n_comp_inner; ++di)
2703 {
2704 Number *out_ = out + di * dofs_per_component;
2705 if (dim > 2)
2706 evaluator.template hessians<2, false, false>(out_, out_);
2707 if (dim > 1)
2708 evaluator.template hessians<1, false, false>(out_, out_);
2709 evaluator.template hessians<0, false, false>(out_, out_);
2710 }
2711
2712 in += dofs_per_component;
2713 out += dofs_per_component;
2714 inv_coefficient += shift_coefficient;
2715 }
2716
2717 return false;
2718 }
2719
2720 private:
2721 template <int n_components>
2722 static inline void
2723 vmult(const Number *inverse_coefficients,
2724 const Number *src,
2725 Number *dst,
2726 const unsigned int dofs_per_component,
2727 const unsigned int n_given_components = 0)
2728 {
2729 const unsigned int n_desired_components =
2730 (n_components > -1) ? n_components : n_given_components;
2731
2732 std::array<Number, dim + 2> tmp = {};
2733 Assert(n_desired_components <= dim + 2,
2734 ExcMessage(
2735 "Number of components larger than dim+2 not supported."));
2736
2737 for (unsigned int d = 0; d < n_desired_components; ++d)
2738 tmp[d] = src[d * dofs_per_component];
2739
2740 for (unsigned int d1 = 0; d1 < n_desired_components; ++d1)
2741 {
2742 const Number *inv_coeff_row =
2743 &inverse_coefficients[d1 * n_desired_components];
2744 Number sum = inv_coeff_row[0] * tmp[0];
2745 for (unsigned int d2 = 1; d2 < n_desired_components; ++d2)
2746 sum += inv_coeff_row[d2] * tmp[d2];
2747 dst[d1 * dofs_per_component] = sum;
2748 }
2749 }
2750 };
2751
2752
2753
2760 template <int dim, typename Number>
2761 struct CellwiseInverseMassMatrixImplTransformFromQPoints
2762 {
2763 template <int fe_degree, int n_q_points_1d>
2764 static bool
2765 run(const unsigned int n_desired_components,
2766 const FEEvaluationData<dim, Number, false> &fe_eval,
2767 const Number *in_array,
2768 Number *out_array)
2769 {
2770 static const bool do_inplace =
2771 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
2772
2773 Assert(fe_eval.get_shape_info().element_type !=
2774 MatrixFreeFunctions::ElementType::tensor_none,
2775 ExcNotImplemented());
2776
2777 const auto &inverse_shape =
2778 do_inplace ?
2779 fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
2780 fe_eval.get_shape_info().data.front().inverse_shape_values;
2781
2782 const std::size_t dofs_per_component =
2783 do_inplace ? Utilities::pow(fe_degree + 1, dim) :
2784 fe_eval.get_shape_info().dofs_per_component_on_cell;
2785 const std::size_t n_q_points = do_inplace ?
2786 Utilities::pow(fe_degree + 1, dim) :
2787 fe_eval.get_shape_info().n_q_points;
2788
2789 using Number2 =
2790 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2791 EvaluatorTensorProduct<do_inplace ? evaluate_evenodd : evaluate_general,
2792 dim,
2793 fe_degree + 1,
2794 n_q_points_1d,
2795 Number,
2796 Number2>
2797 evaluator({},
2798 {},
2799 inverse_shape,
2800 fe_eval.get_shape_info().data.front().fe_degree + 1,
2801 fe_eval.get_shape_info().data.front().n_q_points_1d);
2802
2803 for (unsigned int d = 0; d < n_desired_components; ++d)
2804 {
2805 const Number *in = in_array + d * n_q_points;
2806 Number *out = out_array + d * dofs_per_component;
2807
2808 auto *temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
2809 auto *temp_2 = do_inplace ?
2810 out :
2811 (temp_1 + std::max(n_q_points, dofs_per_component));
2812
2813 if (dim == 3)
2814 {
2815 evaluator.template hessians<2, false, false>(in, temp_1);
2816 evaluator.template hessians<1, false, false>(temp_1, temp_2);
2817 evaluator.template hessians<0, false, false>(temp_2, out);
2818 }
2819 if (dim == 2)
2820 {
2821 evaluator.template hessians<1, false, false>(in, temp_1);
2822 evaluator.template hessians<0, false, false>(temp_1, out);
2823 }
2824 if (dim == 1)
2825 evaluator.template hessians<0, false, false>(in, out);
2826 }
2827 return false;
2828 }
2829 };
2830} // end of namespace internal
2831
2832
2833DEAL_II_NAMESPACE_CLOSE
2834
2835#endif
AlignedVector::size
size_type size() const
ArrayView::begin
iterator begin() const
Definition array_view.h:755
ArrayView::data
value_type * data() const noexcept
Definition array_view.h:714
ArrayView::size
std::size_t size() const
Definition array_view.h:737
FEEvaluationData::shape_info_number_type
ScalarNumber shape_info_number_type
Definition fe_evaluation_data.h:127
FEEvaluationData::get_shape_info
const ShapeInfoType & get_shape_info() const
FEEvaluationData::JxW
Number JxW(const unsigned int q_point) const
FEEvaluationData::begin_gradients
const Number * begin_gradients() const
FEEvaluationData::get_scratch_data
ArrayView< Number > get_scratch_data() const
FEEvaluationData::begin_values
const Number * begin_values() const
FEEvaluationData::begin_hessians
const Number * begin_hessians() const
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:166
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
evaluation_flags.h
evaluation_kernels_common.h
fe_evaluation_data.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
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 exception_macros.h:647
internal::MatrixFreeFunctions::ElementType
ElementType
Definition shape_info.h:56
internal::MatrixFreeFunctions::tensor_symmetric_no_collocation
@ tensor_symmetric_no_collocation
Definition shape_info.h:85
internal::MatrixFreeFunctions::tensor_raviart_thomas
@ tensor_raviart_thomas
Definition shape_info.h:110
internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0
@ tensor_symmetric_plus_dg0
Definition shape_info.h:103
internal::MatrixFreeFunctions::truncated_tensor
@ truncated_tensor
Definition shape_info.h:96
internal::MatrixFreeFunctions::tensor_none
@ tensor_none
Definition shape_info.h:122
internal::MatrixFreeFunctions::tensor_nedelec
@ tensor_nedelec
Definition shape_info.h:117
internal::MatrixFreeFunctions::tensor_general
@ tensor_general
Definition shape_info.h:90
data
std::vector< index_type > data
Definition mpi.cc:734
EvaluationFlags::EvaluationFlags
EvaluationFlags
The EvaluationFlags enum.
Definition evaluation_flags.h:42
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::evaluate_hessians_collocation
void evaluate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:1396
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition tensor_product_kernels.h:2451
internal::integrate_gradients_collocation
void integrate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, Number *values, const Number *gradients, const bool add_into_values_array)
Definition evaluation_kernels.h:1333
internal::evaluate_hessians_slow
void evaluate_hessians_slow(const unsigned int n_components, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:1546
internal::apply_matrix_vector_product
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
Definition tensor_product_kernels.h:116
internal::integrate_hessians_collocation
void integrate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
Definition evaluation_kernels.h:1469
internal::evaluate_gradients_collocation
void evaluate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, const Number *values, Number *gradients)
Definition evaluation_kernels.h:1279
internal::integrate_hessians_slow
void integrate_hessians_slow(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, Number *values_dofs, const bool add_into_values_array)
Definition evaluation_kernels.h:1640
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
shape_info.h
internal::CellwiseInverseMassMatrixImplBasic::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:2523
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)
Definition evaluation_kernels.h:2527
internal::CellwiseInverseMassMatrixImplFlexible::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:2599
internal::CellwiseInverseMassMatrixImplFlexible::run
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const ArrayView< const Number > &inverse_coefficients, const bool dyadic_coefficients, const Number *in_array, Number *out_array)
Definition evaluation_kernels.h:2603
internal::CellwiseInverseMassMatrixImplFlexible::vmult
static void vmult(const Number *inverse_coefficients, const Number *src, Number *dst, const unsigned int dofs_per_component, const unsigned int n_given_components=0)
Definition evaluation_kernels.h:2723
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:2765
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:1047
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:926
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:1197
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:1766
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:1789
internal::FEEvaluationImplCollocation::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:1757
internal::FEEvaluationImplSelector::evaluate_or_integrate
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
Definition evaluation_kernels.h:2497
internal::FEEvaluationImplSelector::evaluate_or_integrate
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array, std::bool_constant< true >)
Definition evaluation_kernels.h:2480
internal::FEEvaluationImplSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array_in=false)
Definition evaluation_kernels.h:2258
internal::FEEvaluationImplSelector::evaluate_or_integrate
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array, std::bool_constant< false >)
Definition evaluation_kernels.h:2465
internal::FEEvaluationImplTransformToCollocation::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:1875
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:1838
internal::FEEvaluationImpl< MatrixFreeFunctions::tensor_nedelec, dim, fe_degree, n_q_points_1d, Number >::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:2084
internal::FEEvaluationImpl< MatrixFreeFunctions::tensor_raviart_thomas, dim, fe_degree, n_q_points_1d, Number >::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:1932
internal::FEEvaluationImpl::variant
static const EvaluatorVariant variant
Definition evaluation_kernels.h:113
internal::FEEvaluationImpl::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:116
internal::FEEvaluationImpl::Eval
EvaluatorTensorProduct< variant, dim, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
Definition evaluation_kernels.h:123
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:373
internal::FEEvaluationImpl::create_evaluator_tensor_product
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > *univariate_shape_data)
Definition evaluation_kernels.h:139
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:193
internal::MatrixFreeFunctions::ShapeInfo::data
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:490
tensor_product_kernels.h
Generated by doxygen 1.9.8