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