Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
evaluation_kernels.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_evaluation_kernels_h
18 #define dealii_matrix_free_evaluation_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/utilities.h>
23 #include <deal.II/base/vectorization.h>
24 
25 #include <deal.II/matrix_free/shape_info.h>
26 #include <deal.II/matrix_free/tensor_product_kernels.h>
27 
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 
32 
33 namespace 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 
104  template <MatrixFreeFunctions::ElementType type,
105  int dim,
106  int fe_degree,
107  int n_q_points_1d,
108  int n_components,
109  typename Number>
110  struct FEEvaluationImpl
111  {
112  static void
113  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
114  const Number * values_dofs_actual,
115  Number * values_quad,
116  Number * gradients_quad,
117  Number * hessians_quad,
118  Number * scratch_data,
119  const bool evaluate_values,
120  const bool evaluate_gradients,
121  const bool evaluate_hessians);
122 
123  static void
124  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
125  Number * values_dofs_actual,
126  Number * values_quad,
127  Number * gradients_quad,
128  Number * scratch_data,
129  const bool integrate_values,
130  const bool integrate_gradients,
131  const bool add_into_values_array);
132  };
133 
134 
135 
136  template <MatrixFreeFunctions::ElementType type,
137  int dim,
138  int fe_degree,
139  int n_q_points_1d,
140  int n_components,
141  typename Number>
142  inline void
143  FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, n_components, Number>::
144  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
145  const Number * values_dofs_actual,
146  Number * values_quad,
147  Number * gradients_quad,
148  Number * hessians_quad,
149  Number * scratch_data,
150  const bool evaluate_values,
151  const bool evaluate_gradients,
152  const bool evaluate_hessians)
153  {
154  if (evaluate_values == false && evaluate_gradients == false &&
155  evaluate_hessians == false)
156  return;
157 
158  const EvaluatorVariant variant =
159  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
160  using Eval = EvaluatorTensorProduct<variant,
161  dim,
162  fe_degree + 1,
163  n_q_points_1d,
164  Number>;
165  Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
166  shape_info.shape_values,
167  variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
168  shape_info.shape_gradients,
169  variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
170  shape_info.shape_hessians,
171  shape_info.fe_degree + 1,
172  shape_info.n_q_points_1d);
173 
174  const unsigned int temp_size =
175  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
176  0 :
177  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
178  Eval::n_rows_of_product :
179  Eval::n_columns_of_product);
180  Number *temp1;
181  Number *temp2;
182  if (temp_size == 0)
183  {
184  temp1 = scratch_data;
185  temp2 = temp1 +
186  std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
187  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
188  }
189  else
190  {
191  temp1 = scratch_data;
192  temp2 = temp1 + temp_size;
193  }
194 
195  const unsigned int n_q_points =
196  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
197  const unsigned int dofs_per_comp =
198  (type == MatrixFreeFunctions::truncated_tensor) ?
199  Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
200  shape_info.dofs_per_component_on_cell;
201  const Number *values_dofs = values_dofs_actual;
202  if (type == MatrixFreeFunctions::truncated_tensor)
203  {
204  Number *values_dofs_tmp =
205  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
206  shape_info.n_q_points));
207  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
208  unsigned int count_p = 0, count_q = 0;
209  for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
210  {
211  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
212  {
213  for (int k = 0; k < degree + 1 - j - i;
214  ++k, ++count_p, ++count_q)
215  for (unsigned int c = 0; c < n_components; ++c)
216  values_dofs_tmp[c * dofs_per_comp + count_q] =
217  values_dofs_actual
218  [c * shape_info.dofs_per_component_on_cell + count_p];
219  for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q)
220  for (unsigned int c = 0; c < n_components; ++c)
221  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
222  }
223  for (int j = degree + 1 - i; j < degree + 1; ++j)
224  for (int k = 0; k < degree + 1; ++k, ++count_q)
225  for (unsigned int c = 0; c < n_components; ++c)
226  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
227  }
228  AssertDimension(count_q, dofs_per_comp);
229  values_dofs = values_dofs_tmp;
230  }
231 
232  switch (dim)
233  {
234  case 1:
235  for (unsigned int c = 0; c < n_components; c++)
236  {
237  if (evaluate_values == true)
238  eval.template values<0, true, false>(values_dofs, values_quad);
239  if (evaluate_gradients == true)
240  eval.template gradients<0, true, false>(values_dofs,
241  gradients_quad);
242  if (evaluate_hessians == true)
243  eval.template hessians<0, true, false>(values_dofs,
244  hessians_quad);
245 
246  // advance the next component in 1D array
247  values_dofs += dofs_per_comp;
248  values_quad += n_q_points;
249  gradients_quad += n_q_points;
250  hessians_quad += n_q_points;
251  }
252  break;
253 
254  case 2:
255  for (unsigned int c = 0; c < n_components; c++)
256  {
257  // grad x
258  if (evaluate_gradients == true)
259  {
260  eval.template gradients<0, true, false>(values_dofs, temp1);
261  eval.template values<1, true, false>(temp1, gradients_quad);
262  }
263  if (evaluate_hessians == true)
264  {
265  // grad xy
266  if (evaluate_gradients == false)
267  eval.template gradients<0, true, false>(values_dofs, temp1);
268  eval.template gradients<1, true, false>(temp1,
269  hessians_quad +
270  2 * n_q_points);
271 
272  // grad xx
273  eval.template hessians<0, true, false>(values_dofs, temp1);
274  eval.template values<1, true, false>(temp1, hessians_quad);
275  }
276 
277  // grad y
278  eval.template values<0, true, false>(values_dofs, temp1);
279  if (evaluate_gradients == true)
280  eval.template gradients<1, true, false>(temp1,
281  gradients_quad +
282  n_q_points);
283 
284  // grad yy
285  if (evaluate_hessians == true)
286  eval.template hessians<1, true, false>(temp1,
287  hessians_quad +
288  n_q_points);
289 
290  // val: can use values applied in x
291  if (evaluate_values == true)
292  eval.template values<1, true, false>(temp1, values_quad);
293 
294  // advance to the next component in 1D array
295  values_dofs += dofs_per_comp;
296  values_quad += n_q_points;
297  gradients_quad += 2 * n_q_points;
298  hessians_quad += 3 * n_q_points;
299  }
300  break;
301 
302  case 3:
303  for (unsigned int c = 0; c < n_components; c++)
304  {
305  if (evaluate_gradients == true)
306  {
307  // grad x
308  eval.template gradients<0, true, false>(values_dofs, temp1);
309  eval.template values<1, true, false>(temp1, temp2);
310  eval.template values<2, true, false>(temp2, gradients_quad);
311  }
312 
313  if (evaluate_hessians == true)
314  {
315  // grad xz
316  if (evaluate_gradients == false)
317  {
318  eval.template gradients<0, true, false>(values_dofs,
319  temp1);
320  eval.template values<1, true, false>(temp1, temp2);
321  }
322  eval.template gradients<2, true, false>(temp2,
323  hessians_quad +
324  4 * n_q_points);
325 
326  // grad xy
327  eval.template gradients<1, true, false>(temp1, temp2);
328  eval.template values<2, true, false>(temp2,
329  hessians_quad +
330  3 * n_q_points);
331 
332  // grad xx
333  eval.template hessians<0, true, false>(values_dofs, temp1);
334  eval.template values<1, true, false>(temp1, temp2);
335  eval.template values<2, true, false>(temp2, hessians_quad);
336  }
337 
338  // grad y
339  eval.template values<0, true, false>(values_dofs, temp1);
340  if (evaluate_gradients == true)
341  {
342  eval.template gradients<1, true, false>(temp1, temp2);
343  eval.template values<2, true, false>(temp2,
344  gradients_quad +
345  n_q_points);
346  }
347 
348  if (evaluate_hessians == true)
349  {
350  // grad yz
351  if (evaluate_gradients == false)
352  eval.template gradients<1, true, false>(temp1, temp2);
353  eval.template gradients<2, true, false>(temp2,
354  hessians_quad +
355  5 * n_q_points);
356 
357  // grad yy
358  eval.template hessians<1, true, false>(temp1, temp2);
359  eval.template values<2, true, false>(temp2,
360  hessians_quad +
361  n_q_points);
362  }
363 
364  // grad z: can use the values applied in x direction stored in
365  // temp1
366  eval.template values<1, true, false>(temp1, temp2);
367  if (evaluate_gradients == true)
368  eval.template gradients<2, true, false>(temp2,
369  gradients_quad +
370  2 * n_q_points);
371 
372  // grad zz: can use the values applied in x and y direction stored
373  // in temp2
374  if (evaluate_hessians == true)
375  eval.template hessians<2, true, false>(temp2,
376  hessians_quad +
377  2 * n_q_points);
378 
379  // val: can use the values applied in x & y direction stored in
380  // temp2
381  if (evaluate_values == true)
382  eval.template values<2, true, false>(temp2, values_quad);
383 
384  // advance to the next component in 1D array
385  values_dofs += dofs_per_comp;
386  values_quad += n_q_points;
387  gradients_quad += 3 * n_q_points;
388  hessians_quad += 6 * n_q_points;
389  }
390  break;
391 
392  default:
393  AssertThrow(false, ExcNotImplemented());
394  }
395 
396  // case additional dof for FE_Q_DG0: add values; gradients and second
397  // derivatives evaluate to zero
398  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 &&
399  evaluate_values)
400  {
401  values_quad -= n_components * n_q_points;
402  values_dofs -= n_components * dofs_per_comp;
403  for (unsigned int c = 0; c < n_components; ++c)
404  for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
405  values_quad[c * shape_info.n_q_points + q] +=
406  values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
407  }
408  }
409 
410 
411 
412  template <MatrixFreeFunctions::ElementType type,
413  int dim,
414  int fe_degree,
415  int n_q_points_1d,
416  int n_components,
417  typename Number>
418  inline void
419  FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, n_components, Number>::
420  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
421  Number * values_dofs_actual,
422  Number * values_quad,
423  Number * gradients_quad,
424  Number * scratch_data,
425  const bool integrate_values,
426  const bool integrate_gradients,
427  const bool add_into_values_array)
428  {
429  const EvaluatorVariant variant =
430  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
431  using Eval = EvaluatorTensorProduct<variant,
432  dim,
433  fe_degree + 1,
434  n_q_points_1d,
435  Number>;
436  Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
437  shape_info.shape_values,
438  variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
439  shape_info.shape_gradients,
440  variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
441  shape_info.shape_hessians,
442  shape_info.fe_degree + 1,
443  shape_info.n_q_points_1d);
444 
445  const unsigned int temp_size =
446  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
447  0 :
448  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
449  Eval::n_rows_of_product :
450  Eval::n_columns_of_product);
451  Number *temp1;
452  Number *temp2;
453  if (temp_size == 0)
454  {
455  temp1 = scratch_data;
456  temp2 = temp1 +
457  std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
458  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
459  }
460  else
461  {
462  temp1 = scratch_data;
463  temp2 = temp1 + temp_size;
464  }
465 
466  const unsigned int n_q_points =
467  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
468  const unsigned int dofs_per_comp =
469  (type == MatrixFreeFunctions::truncated_tensor) ?
470  Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
471  shape_info.dofs_per_component_on_cell;
472  // expand dof_values to tensor product for truncated tensor products
473  Number *values_dofs =
474  (type == MatrixFreeFunctions::truncated_tensor) ?
475  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
476  shape_info.n_q_points)) :
477  values_dofs_actual;
478 
479  switch (dim)
480  {
481  case 1:
482  for (unsigned int c = 0; c < n_components; c++)
483  {
484  if (integrate_values == true)
485  {
486  if (add_into_values_array == false)
487  eval.template values<0, false, false>(values_quad,
488  values_dofs);
489  else
490  eval.template values<0, false, true>(values_quad,
491  values_dofs);
492  }
493  if (integrate_gradients == true)
494  {
495  if (integrate_values == true || add_into_values_array == true)
496  eval.template gradients<0, false, true>(gradients_quad,
497  values_dofs);
498  else
499  eval.template gradients<0, false, false>(gradients_quad,
500  values_dofs);
501  }
502 
503  // advance to the next component in 1D array
504  values_dofs += dofs_per_comp;
505  values_quad += n_q_points;
506  gradients_quad += n_q_points;
507  }
508  break;
509 
510  case 2:
511  for (unsigned int c = 0; c < n_components; c++)
512  {
513  if (integrate_values == true && integrate_gradients == false)
514  {
515  eval.template values<1, false, false>(values_quad, temp1);
516  if (add_into_values_array == false)
517  eval.template values<0, false, false>(temp1, values_dofs);
518  else
519  eval.template values<0, false, true>(temp1, values_dofs);
520  }
521  if (integrate_gradients == true)
522  {
523  eval.template gradients<1, false, false>(gradients_quad +
524  n_q_points,
525  temp1);
526  if (integrate_values)
527  eval.template values<1, false, true>(values_quad, temp1);
528  if (add_into_values_array == false)
529  eval.template values<0, false, false>(temp1, values_dofs);
530  else
531  eval.template values<0, false, true>(temp1, values_dofs);
532  eval.template values<1, false, false>(gradients_quad, temp1);
533  eval.template gradients<0, false, true>(temp1, values_dofs);
534  }
535 
536  // advance to the next component in 1D array
537  values_dofs += dofs_per_comp;
538  values_quad += n_q_points;
539  gradients_quad += 2 * n_q_points;
540  }
541  break;
542 
543  case 3:
544  for (unsigned int c = 0; c < n_components; c++)
545  {
546  if (integrate_values == true && integrate_gradients == false)
547  {
548  eval.template values<2, false, false>(values_quad, temp1);
549  eval.template values<1, false, false>(temp1, temp2);
550  if (add_into_values_array == false)
551  eval.template values<0, false, false>(temp2, values_dofs);
552  else
553  eval.template values<0, false, true>(temp2, values_dofs);
554  }
555  if (integrate_gradients == true)
556  {
557  eval.template gradients<2, false, false>(gradients_quad +
558  2 * n_q_points,
559  temp1);
560  if (integrate_values)
561  eval.template values<2, false, true>(values_quad, temp1);
562  eval.template values<1, false, false>(temp1, temp2);
563  eval.template values<2, false, false>(gradients_quad +
564  n_q_points,
565  temp1);
566  eval.template gradients<1, false, true>(temp1, temp2);
567  if (add_into_values_array == false)
568  eval.template values<0, false, false>(temp2, values_dofs);
569  else
570  eval.template values<0, false, true>(temp2, values_dofs);
571  eval.template values<2, false, false>(gradients_quad, temp1);
572  eval.template values<1, false, false>(temp1, temp2);
573  eval.template gradients<0, false, true>(temp2, values_dofs);
574  }
575 
576  // advance to the next component in 1D array
577  values_dofs += dofs_per_comp;
578  values_quad += n_q_points;
579  gradients_quad += 3 * n_q_points;
580  }
581  break;
582 
583  default:
584  AssertThrow(false, ExcNotImplemented());
585  }
586 
587  // case FE_Q_DG0: add values, gradients and second derivatives are zero
588  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
589  {
590  values_dofs -= n_components * dofs_per_comp -
591  shape_info.dofs_per_component_on_cell + 1;
592  values_quad -= n_components * n_q_points;
593  if (integrate_values)
594  for (unsigned int c = 0; c < n_components; ++c)
595  {
596  values_dofs[0] = values_quad[0];
597  for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
598  values_dofs[0] += values_quad[q];
599  values_dofs += dofs_per_comp;
600  values_quad += n_q_points;
601  }
602  else
603  {
604  for (unsigned int c = 0; c < n_components; ++c)
605  values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
606  values_dofs += n_components * shape_info.dofs_per_component_on_cell;
607  }
608  }
609 
610  if (type == MatrixFreeFunctions::truncated_tensor)
611  {
612  values_dofs -= dofs_per_comp * n_components;
613  unsigned int count_p = 0, count_q = 0;
614  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
615  for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
616  {
617  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
618  {
619  for (int k = 0; k < degree + 1 - j - i;
620  ++k, ++count_p, ++count_q)
621  {
622  for (unsigned int c = 0; c < n_components; ++c)
623  values_dofs_actual
624  [c * shape_info.dofs_per_component_on_cell + count_p] =
625  values_dofs[c * dofs_per_comp + count_q];
626  }
627  count_q += j + i;
628  }
629  count_q += i * (degree + 1);
630  }
631  AssertDimension(count_q,
632  Utilities::fixed_power<dim>(shape_info.fe_degree + 1));
633  }
634  }
635 
636 
637 
649  template <EvaluatorVariant variant,
650  int dim,
651  int basis_size_1,
652  int basis_size_2,
653  int n_components,
654  typename Number,
655  typename Number2>
656  struct FEEvaluationImplBasisChange
657  {
658  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
659  "The second dimension must not be smaller than the first");
660 
682 #ifndef DEBUG
683  DEAL_II_ALWAYS_INLINE
684 #endif
685  static void
686  do_forward(
687  const AlignedVector<Number2> &transformation_matrix,
688  const Number * values_in,
689  Number * values_out,
690  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
691  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
692  {
693  Assert(
694  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
695  ExcMessage("The second dimension must not be smaller than the first"));
696 
697  // we do recursion until dim==1 or dim==2 and we have
698  // basis_size_1==basis_size_2. The latter optimization increases
699  // optimization possibilities for the compiler but does only work for
700  // aliased pointers if the sizes are equal.
701  constexpr int next_dim =
702  (dim > 2 ||
703  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
704  dim - 1 :
705  dim;
706 
707  EvaluatorTensorProduct<variant,
708  dim,
709  basis_size_1,
710  (basis_size_1 == 0 ? 0 : basis_size_2),
711  Number,
712  Number2>
713  eval_val(transformation_matrix,
714  AlignedVector<Number2>(),
715  AlignedVector<Number2>(),
716  basis_size_1_variable,
717  basis_size_2_variable);
718  const unsigned int np_1 =
719  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
720  const unsigned int np_2 =
721  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
722  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
723  ExcMessage("Cannot transform with 0-point basis"));
724  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
725  ExcMessage("Cannot transform with 0-point basis"));
726 
727  // run loop backwards to ensure correctness if values_in aliases with
728  // values_out in case with basis_size_1 < basis_size_2
729  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
730  values_out =
731  values_out + n_components * Utilities::fixed_power<dim>(np_2);
732  for (unsigned int c = n_components; c != 0; --c)
733  {
734  values_in -= Utilities::fixed_power<dim>(np_1);
735  values_out -= Utilities::fixed_power<dim>(np_2);
736  if (next_dim < dim)
737  for (unsigned int q = np_1; q != 0; --q)
738  FEEvaluationImplBasisChange<
739  variant,
740  next_dim,
741  basis_size_1,
742  basis_size_2,
743  1,
744  Number,
745  Number2>::do_forward(transformation_matrix,
746  values_in +
747  (q - 1) *
748  Utilities::fixed_power<next_dim>(np_1),
749  values_out +
750  (q - 1) *
751  Utilities::fixed_power<next_dim>(np_2),
752  basis_size_1_variable,
753  basis_size_2_variable);
754 
755  // the recursion stops if dim==1 or if dim==2 and
756  // basis_size_1==basis_size_2 (the latter is used because the
757  // compiler generates nicer code)
758  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
759  {
760  eval_val.template values<0, true, false>(values_in, values_out);
761  eval_val.template values<1, true, false>(values_out, values_out);
762  }
763  else if (dim == 1)
764  eval_val.template values<dim - 1, true, false>(values_in,
765  values_out);
766  else
767  eval_val.template values<dim - 1, true, false>(values_out,
768  values_out);
769  }
770  }
771 
801 #ifndef DEBUG
802  DEAL_II_ALWAYS_INLINE
803 #endif
804  static void
805  do_backward(
806  const AlignedVector<Number2> &transformation_matrix,
807  const bool add_into_result,
808  Number * values_in,
809  Number * values_out,
810  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
811  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
812  {
813  Assert(
814  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
815  ExcMessage("The second dimension must not be smaller than the first"));
816  Assert(add_into_result == false || values_in != values_out,
817  ExcMessage(
818  "Input and output cannot alias with each other when "
819  "adding the result of the basis change to existing data"));
820 
821  constexpr int next_dim =
822  (dim > 2 ||
823  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
824  dim - 1 :
825  dim;
826  EvaluatorTensorProduct<variant,
827  dim,
828  basis_size_1,
829  (basis_size_1 == 0 ? 0 : basis_size_2),
830  Number,
831  Number2>
832  eval_val(transformation_matrix,
833  AlignedVector<Number2>(),
834  AlignedVector<Number2>(),
835  basis_size_1_variable,
836  basis_size_2_variable);
837  const unsigned int np_1 =
838  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
839  const unsigned int np_2 =
840  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
841  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
842  ExcMessage("Cannot transform with 0-point basis"));
843  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
844  ExcMessage("Cannot transform with 0-point basis"));
845 
846  for (unsigned int c = 0; c < n_components; ++c)
847  {
848  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
849  {
850  eval_val.template values<1, false, false>(values_in, values_in);
851  if (add_into_result)
852  eval_val.template values<0, false, true>(values_in, values_out);
853  else
854  eval_val.template values<0, false, false>(values_in,
855  values_out);
856  }
857  else
858  {
859  if (dim == 1 && add_into_result)
860  eval_val.template values<0, false, true>(values_in, values_out);
861  else if (dim == 1)
862  eval_val.template values<0, false, false>(values_in,
863  values_out);
864  else
865  eval_val.template values<dim - 1, false, false>(values_in,
866  values_in);
867  }
868  if (next_dim < dim)
869  for (unsigned int q = 0; q < np_1; ++q)
870  FEEvaluationImplBasisChange<variant,
871  next_dim,
872  basis_size_1,
873  basis_size_2,
874  1,
875  Number,
876  Number2>::
877  do_backward(transformation_matrix,
878  add_into_result,
879  values_in +
880  q * Utilities::fixed_power<next_dim>(np_2),
881  values_out +
882  q * Utilities::fixed_power<next_dim>(np_1),
883  basis_size_1_variable,
884  basis_size_2_variable);
885 
886  values_in += Utilities::fixed_power<dim>(np_2);
887  values_out += Utilities::fixed_power<dim>(np_1);
888  }
889  }
890 
910  static void
911  do_mass(const AlignedVector<Number2> &transformation_matrix,
912  const AlignedVector<Number> & coefficients,
913  const Number * values_in,
914  Number * scratch_data,
915  Number * values_out)
916  {
917  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
918  Number * my_scratch =
919  basis_size_1 != basis_size_2 ? scratch_data : values_out;
920 
921  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
922  Assert(coefficients.size() == size_per_component ||
923  coefficients.size() == n_components * size_per_component,
924  ExcDimensionMismatch(coefficients.size(), size_per_component));
925  const unsigned int stride =
926  coefficients.size() == size_per_component ? 0 : 1;
927 
928  for (unsigned int q = basis_size_1; q != 0; --q)
929  FEEvaluationImplBasisChange<
930  variant,
931  next_dim,
932  basis_size_1,
933  basis_size_2,
934  n_components,
935  Number,
936  Number2>::do_forward(transformation_matrix,
937  values_in +
938  (q - 1) *
939  Utilities::pow(basis_size_1, dim - 1),
940  my_scratch +
941  (q - 1) *
942  Utilities::pow(basis_size_2, dim - 1));
943  EvaluatorTensorProduct<variant,
944  dim,
945  basis_size_1,
946  basis_size_2,
947  Number,
948  Number2>
949  eval_val(transformation_matrix);
950  const unsigned int n_inner_blocks =
951  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
952  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
953  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
954  for (unsigned int c = 0; c < n_components; ++c)
955  {
956  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
957  eval_val.template values_one_line<dim - 1, true, false>(
958  my_scratch + i, my_scratch + i);
959  for (unsigned int q = 0; q < basis_size_2; ++q)
960  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
961  my_scratch[i + q * n_blocks + c * size_per_component] *=
962  coefficients[i + q * n_blocks +
963  c * stride * size_per_component];
964  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
965  eval_val.template values_one_line<dim - 1, false, false>(
966  my_scratch + i, my_scratch + i);
967  }
968  for (unsigned int q = 0; q < basis_size_1; ++q)
969  FEEvaluationImplBasisChange<
970  variant,
971  next_dim,
972  basis_size_1,
973  basis_size_2,
974  n_components,
975  Number,
976  Number2>::do_backward(transformation_matrix,
977  false,
978  my_scratch +
979  q * Utilities::pow(basis_size_2, dim - 1),
980  values_out +
981  q * Utilities::pow(basis_size_1, dim - 1));
982  }
983  };
984 
985 
986 
1001  template <int dim, int fe_degree, int n_components, typename Number>
1002  struct FEEvaluationImplCollocation
1003  {
1004  static void
1005  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1006  const Number * values_dofs,
1007  Number * values_quad,
1008  Number * gradients_quad,
1009  Number * hessians_quad,
1010  Number * scratch_data,
1011  const bool evaluate_values,
1012  const bool evaluate_gradients,
1013  const bool evaluate_hessians);
1014 
1015  static void
1016  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1017  Number * values_dofs,
1018  Number * values_quad,
1019  Number * gradients_quad,
1020  Number * scratch_data,
1021  const bool integrate_values,
1022  const bool integrate_gradients,
1023  const bool add_into_values_array);
1024  };
1025 
1026 
1027 
1028  template <int dim, int fe_degree, int n_components, typename Number>
1029  inline void
1030  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::evaluate(
1031  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1032  const Number * values_dofs,
1033  Number * values_quad,
1034  Number * gradients_quad,
1035  Number * hessians_quad,
1036  Number *,
1037  const bool evaluate_values,
1038  const bool evaluate_gradients,
1039  const bool evaluate_hessians)
1040  {
1041  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1042  (fe_degree + 2) / 2 * (fe_degree + 1));
1043 
1044  EvaluatorTensorProduct<evaluate_evenodd,
1045  dim,
1046  fe_degree + 1,
1047  fe_degree + 1,
1048  Number>
1049  eval(AlignedVector<Number>(),
1050  shape_info.shape_gradients_collocation_eo,
1051  shape_info.shape_hessians_collocation_eo);
1052  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1053 
1054  for (unsigned int c = 0; c < n_components; c++)
1055  {
1056  if (evaluate_values == true)
1057  for (unsigned int i = 0; i < n_q_points; ++i)
1058  values_quad[i] = values_dofs[i];
1059  if (evaluate_gradients == true || evaluate_hessians == true)
1060  {
1061  eval.template gradients<0, true, false>(values_dofs,
1062  gradients_quad);
1063  if (dim > 1)
1064  eval.template gradients<1, true, false>(values_dofs,
1065  gradients_quad +
1066  n_q_points);
1067  if (dim > 2)
1068  eval.template gradients<2, true, false>(values_dofs,
1069  gradients_quad +
1070  2 * n_q_points);
1071  }
1072  if (evaluate_hessians == true)
1073  {
1074  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1075  if (dim > 1)
1076  {
1077  eval.template gradients<1, true, false>(gradients_quad,
1078  hessians_quad +
1079  dim * n_q_points);
1080  eval.template hessians<1, true, false>(values_dofs,
1081  hessians_quad +
1082  n_q_points);
1083  }
1084  if (dim > 2)
1085  {
1086  eval.template gradients<2, true, false>(gradients_quad,
1087  hessians_quad +
1088  4 * n_q_points);
1089  eval.template gradients<2, true, false>(
1090  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
1091  eval.template hessians<2, true, false>(values_dofs,
1092  hessians_quad +
1093  2 * n_q_points);
1094  }
1095  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1096  }
1097  gradients_quad += dim * n_q_points;
1098  values_quad += n_q_points;
1099  values_dofs += n_q_points;
1100  }
1101  }
1102 
1103 
1104 
1105  template <int dim, int fe_degree, int n_components, typename Number>
1106  inline void
1107  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::integrate(
1108  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1109  Number * values_dofs,
1110  Number * values_quad,
1111  Number * gradients_quad,
1112  Number *,
1113  const bool integrate_values,
1114  const bool integrate_gradients,
1115  const bool add_into_values_array)
1116  {
1117  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1118  (fe_degree + 2) / 2 * (fe_degree + 1));
1119 
1120  EvaluatorTensorProduct<evaluate_evenodd,
1121  dim,
1122  fe_degree + 1,
1123  fe_degree + 1,
1124  Number>
1125  eval(AlignedVector<Number>(),
1126  shape_info.shape_gradients_collocation_eo,
1127  shape_info.shape_hessians_collocation_eo);
1128  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1129 
1130  for (unsigned int c = 0; c < n_components; c++)
1131  {
1132  if (integrate_values == true && add_into_values_array == false)
1133  for (unsigned int i = 0; i < n_q_points; ++i)
1134  values_dofs[i] = values_quad[i];
1135  else if (integrate_values == true)
1136  for (unsigned int i = 0; i < n_q_points; ++i)
1137  values_dofs[i] += values_quad[i];
1138  if (integrate_gradients == true)
1139  {
1140  if (integrate_values == true || add_into_values_array == true)
1141  eval.template gradients<0, false, true>(gradients_quad,
1142  values_dofs);
1143  else
1144  eval.template gradients<0, false, false>(gradients_quad,
1145  values_dofs);
1146  if (dim > 1)
1147  eval.template gradients<1, false, true>(gradients_quad +
1148  n_q_points,
1149  values_dofs);
1150  if (dim > 2)
1151  eval.template gradients<2, false, true>(gradients_quad +
1152  2 * n_q_points,
1153  values_dofs);
1154  }
1155  gradients_quad += dim * n_q_points;
1156  values_quad += n_q_points;
1157  values_dofs += n_q_points;
1158  }
1159  }
1160 
1161 
1162 
1175  template <int dim,
1176  int fe_degree,
1177  int n_q_points_1d,
1178  int n_components,
1179  typename Number>
1180  struct FEEvaluationImplTransformToCollocation
1181  {
1182  static void
1183  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1184  const Number * values_dofs,
1185  Number * values_quad,
1186  Number * gradients_quad,
1187  Number * hessians_quad,
1188  Number * scratch_data,
1189  const bool evaluate_values,
1190  const bool evaluate_gradients,
1191  const bool evaluate_hessians);
1192 
1193  static void
1194  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1195  Number * values_dofs,
1196  Number * values_quad,
1197  Number * gradients_quad,
1198  Number * scratch_data,
1199  const bool integrate_values,
1200  const bool integrate_gradients,
1201  const bool add_into_values_array);
1202  };
1203 
1204 
1205 
1206  template <int dim,
1207  int fe_degree,
1208  int n_q_points_1d,
1209  int n_components,
1210  typename Number>
1211  inline void
1212  FEEvaluationImplTransformToCollocation<
1213  dim,
1214  fe_degree,
1215  n_q_points_1d,
1216  n_components,
1217  Number>::evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1218  const Number * values_dofs,
1219  Number * values_quad,
1220  Number *gradients_quad,
1221  Number *hessians_quad,
1222  Number *,
1223  const bool,
1224  const bool evaluate_gradients,
1225  const bool evaluate_hessians)
1226  {
1227  Assert(n_q_points_1d > fe_degree,
1228  ExcMessage("You lose information when going to a collocation space "
1229  "of lower degree, so the evaluation results would be "
1230  "wrong. Thus, this class does not permit the desired "
1231  "operation."));
1232  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1233 
1234  for (unsigned int c = 0; c < n_components; c++)
1235  {
1236  FEEvaluationImplBasisChange<
1237  evaluate_evenodd,
1238  dim,
1239  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1240  n_q_points_1d,
1241  1,
1242  Number,
1243  Number>::do_forward(shape_info.shape_values_eo,
1244  values_dofs,
1245  values_quad);
1246 
1247  // apply derivatives in the collocation space
1248  if (evaluate_gradients == true || evaluate_hessians == true)
1249  FEEvaluationImplCollocation<dim, n_q_points_1d - 1, 1, Number>::
1250  evaluate(shape_info,
1251  values_quad,
1252  nullptr,
1253  gradients_quad,
1254  hessians_quad,
1255  nullptr,
1256  false,
1257  evaluate_gradients,
1258  evaluate_hessians);
1259 
1260  values_dofs += shape_info.dofs_per_component_on_cell;
1261  values_quad += n_q_points;
1262  gradients_quad += dim * n_q_points;
1263  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1264  }
1265  }
1266 
1267 
1268 
1269  template <int dim,
1270  int fe_degree,
1271  int n_q_points_1d,
1272  int n_components,
1273  typename Number>
1274  inline void
1275  FEEvaluationImplTransformToCollocation<
1276  dim,
1277  fe_degree,
1278  n_q_points_1d,
1279  n_components,
1280  Number>::integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1281  Number *values_dofs,
1282  Number *values_quad,
1283  Number *gradients_quad,
1284  Number *,
1285  const bool integrate_values,
1286  const bool integrate_gradients,
1287  const bool add_into_values_array)
1288  {
1289  Assert(n_q_points_1d > fe_degree,
1290  ExcMessage("You lose information when going to a collocation space "
1291  "of lower degree, so the evaluation results would be "
1292  "wrong. Thus, this class does not permit the desired "
1293  "operation."));
1294  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1295  (n_q_points_1d + 1) / 2 * n_q_points_1d);
1296  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1297 
1298  for (unsigned int c = 0; c < n_components; c++)
1299  {
1300  // apply derivatives in collocation space
1301  if (integrate_gradients == true)
1302  FEEvaluationImplCollocation<dim, n_q_points_1d - 1, 1, Number>::
1303  integrate(shape_info,
1304  values_quad,
1305  nullptr,
1306  gradients_quad,
1307  nullptr,
1308  false,
1309  integrate_gradients,
1310  /*add_into_values_array=*/integrate_values);
1311 
1312  // transform back to the original space
1313  FEEvaluationImplBasisChange<
1314  evaluate_evenodd,
1315  dim,
1316  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1317  n_q_points_1d,
1318  1,
1319  Number,
1320  Number>::do_backward(shape_info.shape_values_eo,
1321  add_into_values_array,
1322  values_quad,
1323  values_dofs);
1324  gradients_quad += dim * n_q_points;
1325  values_quad += n_q_points;
1326  values_dofs += shape_info.dofs_per_component_on_cell;
1327  }
1328  }
1329 
1330 
1331 
1332  template <bool symmetric_evaluate,
1333  int dim,
1334  int fe_degree,
1335  int n_q_points_1d,
1336  int n_components,
1337  typename Number>
1338  struct FEFaceEvaluationImpl
1339  {
1340  // We enable a transformation to collocation for derivatives if it gives
1341  // correct results (first two conditions), if it is the most efficient
1342  // choice in terms of operation counts (third condition) and if we were
1343  // able to initialize the fields in shape_info.templates.h from the
1344  // polynomials (fourth condition).
1345  static constexpr bool use_collocation =
1346  symmetric_evaluate &&
1347  n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1348  n_q_points_1d < 200;
1349 
1350  static void
1351  evaluate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1352  Number * values_dofs,
1353  Number * values_quad,
1354  Number * gradients_quad,
1355  Number * scratch_data,
1356  const bool evaluate_val,
1357  const bool evaluate_grad,
1358  const unsigned int subface_index)
1359  {
1360  const AlignedVector<Number> &val1 =
1361  symmetric_evaluate ?
1362  data.shape_values_eo :
1363  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1364  data.shape_values :
1365  data.values_within_subface[subface_index % 2]);
1366  const AlignedVector<Number> &val2 =
1367  symmetric_evaluate ?
1368  data.shape_values_eo :
1369  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1370  data.shape_values :
1371  data.values_within_subface[subface_index / 2]);
1372 
1373  const AlignedVector<Number> &grad1 =
1374  symmetric_evaluate ?
1375  data.shape_gradients_eo :
1376  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1377  data.shape_gradients :
1378  data.gradients_within_subface[subface_index % 2]);
1379  const AlignedVector<Number> &grad2 =
1380  symmetric_evaluate ?
1381  data.shape_gradients_eo :
1382  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1383  data.shape_gradients :
1384  data.gradients_within_subface[subface_index / 2]);
1385 
1386  using Eval =
1387  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1388  internal::evaluate_evenodd :
1389  internal::evaluate_general,
1390  dim - 1,
1391  fe_degree + 1,
1392  n_q_points_1d,
1393  Number>;
1394  Eval eval1(val1,
1395  grad1,
1396  AlignedVector<Number>(),
1397  data.fe_degree + 1,
1398  data.n_q_points_1d);
1399  Eval eval2(val2,
1400  grad2,
1401  AlignedVector<Number>(),
1402  data.fe_degree + 1,
1403  data.n_q_points_1d);
1404 
1405  const unsigned int size_deg =
1406  fe_degree > -1 ?
1407  Utilities::pow(fe_degree + 1, dim - 1) :
1408  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
1409 
1410  const unsigned int n_q_points = fe_degree > -1 ?
1411  Utilities::pow(n_q_points_1d, dim - 1) :
1412  data.n_q_points_face;
1413 
1414  if (evaluate_grad == false)
1415  for (unsigned int c = 0; c < n_components; ++c)
1416  {
1417  switch (dim)
1418  {
1419  case 3:
1420  eval1.template values<0, true, false>(values_dofs,
1421  values_quad);
1422  eval2.template values<1, true, false>(values_quad,
1423  values_quad);
1424  break;
1425  case 2:
1426  eval1.template values<0, true, false>(values_dofs,
1427  values_quad);
1428  break;
1429  case 1:
1430  values_quad[0] = values_dofs[0];
1431  break;
1432  default:
1433  Assert(false, ExcNotImplemented());
1434  }
1435  values_dofs += 2 * size_deg;
1436  values_quad += n_q_points;
1437  }
1438  else
1439  for (unsigned int c = 0; c < n_components; ++c)
1440  {
1441  switch (dim)
1442  {
1443  case 3:
1444  if (use_collocation)
1445  {
1446  eval1.template values<0, true, false>(values_dofs,
1447  values_quad);
1448  eval1.template values<1, true, false>(values_quad,
1449  values_quad);
1450  internal::EvaluatorTensorProduct<
1451  internal::evaluate_evenodd,
1452  dim - 1,
1453  n_q_points_1d,
1454  n_q_points_1d,
1455  Number>
1456  eval_grad(AlignedVector<Number>(),
1457  data.shape_gradients_collocation_eo,
1458  AlignedVector<Number>());
1459  eval_grad.template gradients<0, true, false>(
1460  values_quad, gradients_quad);
1461  eval_grad.template gradients<1, true, false>(
1462  values_quad, gradients_quad + n_q_points);
1463  }
1464  else
1465  {
1466  eval1.template gradients<0, true, false>(values_dofs,
1467  scratch_data);
1468  eval2.template values<1, true, false>(scratch_data,
1469  gradients_quad);
1470 
1471  eval1.template values<0, true, false>(values_dofs,
1472  scratch_data);
1473  eval2.template gradients<1, true, false>(scratch_data,
1474  gradients_quad +
1475  n_q_points);
1476 
1477  if (evaluate_val == true)
1478  eval2.template values<1, true, false>(scratch_data,
1479  values_quad);
1480  }
1481  eval1.template values<0, true, false>(values_dofs + size_deg,
1482  scratch_data);
1483  eval2.template values<1, true, false>(
1484  scratch_data, gradients_quad + (dim - 1) * n_q_points);
1485 
1486  break;
1487  case 2:
1488  eval1.template values<0, true, false>(values_dofs + size_deg,
1489  gradients_quad +
1490  (dim - 1) *
1491  n_q_points);
1492  eval1.template gradients<0, true, false>(values_dofs,
1493  gradients_quad);
1494  if (evaluate_val == true)
1495  eval1.template values<0, true, false>(values_dofs,
1496  values_quad);
1497  break;
1498  case 1:
1499  values_quad[0] = values_dofs[0];
1500  gradients_quad[0] = values_dofs[1];
1501  break;
1502  default:
1503  AssertThrow(false, ExcNotImplemented());
1504  }
1505  values_dofs += 2 * size_deg;
1506  values_quad += n_q_points;
1507  gradients_quad += dim * n_q_points;
1508  }
1509  }
1510 
1511  static void
1512  integrate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1513  Number * values_dofs,
1514  Number * values_quad,
1515  Number * gradients_quad,
1516  Number * scratch_data,
1517  const bool integrate_val,
1518  const bool integrate_grad,
1519  const unsigned int subface_index)
1520  {
1521  const AlignedVector<Number> &val1 =
1522  symmetric_evaluate ?
1523  data.shape_values_eo :
1524  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1525  data.shape_values :
1526  data.values_within_subface[subface_index % 2]);
1527  const AlignedVector<Number> &val2 =
1528  symmetric_evaluate ?
1529  data.shape_values_eo :
1530  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1531  data.shape_values :
1532  data.values_within_subface[subface_index / 2]);
1533 
1534  const AlignedVector<Number> &grad1 =
1535  symmetric_evaluate ?
1536  data.shape_gradients_eo :
1537  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1538  data.shape_gradients :
1539  data.gradients_within_subface[subface_index % 2]);
1540  const AlignedVector<Number> &grad2 =
1541  symmetric_evaluate ?
1542  data.shape_gradients_eo :
1543  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1544  data.shape_gradients :
1545  data.gradients_within_subface[subface_index / 2]);
1546 
1547  using Eval =
1548  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1549  internal::evaluate_evenodd :
1550  internal::evaluate_general,
1551  dim - 1,
1552  fe_degree + 1,
1553  n_q_points_1d,
1554  Number>;
1555  Eval eval1(val1, grad1, val1, data.fe_degree + 1, data.n_q_points_1d);
1556  Eval eval2(val2, grad2, val1, data.fe_degree + 1, data.n_q_points_1d);
1557 
1558  const unsigned int size_deg =
1559  fe_degree > -1 ?
1560  Utilities::pow(fe_degree + 1, dim - 1) :
1561  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
1562 
1563  const unsigned int n_q_points = fe_degree > -1 ?
1564  Utilities::pow(n_q_points_1d, dim - 1) :
1565  data.n_q_points_face;
1566 
1567  if (integrate_grad == false)
1568  for (unsigned int c = 0; c < n_components; ++c)
1569  {
1570  switch (dim)
1571  {
1572  case 3:
1573  eval2.template values<1, false, false>(values_quad,
1574  values_quad);
1575  eval1.template values<0, false, false>(values_quad,
1576  values_dofs);
1577  break;
1578  case 2:
1579  eval1.template values<0, false, false>(values_quad,
1580  values_dofs);
1581  break;
1582  case 1:
1583  values_dofs[0] = values_quad[0];
1584  break;
1585  default:
1586  Assert(false, ExcNotImplemented());
1587  }
1588  values_dofs += 2 * size_deg;
1589  values_quad += n_q_points;
1590  }
1591  else
1592  for (unsigned int c = 0; c < n_components; ++c)
1593  {
1594  switch (dim)
1595  {
1596  case 3:
1597  eval2.template values<1, false, false>(gradients_quad +
1598  2 * n_q_points,
1599  gradients_quad +
1600  2 * n_q_points);
1601  eval1.template values<0, false, false>(
1602  gradients_quad + 2 * n_q_points, values_dofs + size_deg);
1603  if (use_collocation)
1604  {
1605  internal::EvaluatorTensorProduct<
1606  internal::evaluate_evenodd,
1607  dim - 1,
1608  n_q_points_1d,
1609  n_q_points_1d,
1610  Number>
1611  eval_grad(AlignedVector<Number>(),
1612  data.shape_gradients_collocation_eo,
1613  AlignedVector<Number>());
1614  if (integrate_val)
1615  eval_grad.template gradients<1, false, true>(
1616  gradients_quad + n_q_points, values_quad);
1617  else
1618  eval_grad.template gradients<1, false, false>(
1619  gradients_quad + n_q_points, values_quad);
1620  eval_grad.template gradients<0, false, true>(
1621  gradients_quad, values_quad);
1622  eval1.template values<1, false, false>(values_quad,
1623  values_quad);
1624  eval1.template values<0, false, false>(values_quad,
1625  values_dofs);
1626  }
1627  else
1628  {
1629  if (integrate_val)
1630  {
1631  eval2.template values<1, false, false>(values_quad,
1632  scratch_data);
1633  eval2.template gradients<1, false, true>(
1634  gradients_quad + n_q_points, scratch_data);
1635  }
1636  else
1637  eval2.template gradients<1, false, false>(
1638  gradients_quad + n_q_points, scratch_data);
1639 
1640  eval1.template values<0, false, false>(scratch_data,
1641  values_dofs);
1642  eval2.template values<1, false, false>(gradients_quad,
1643  scratch_data);
1644  eval1.template gradients<0, false, true>(scratch_data,
1645  values_dofs);
1646  }
1647  break;
1648  case 2:
1649  eval1.template values<0, false, false>(
1650  gradients_quad + n_q_points, values_dofs + size_deg);
1651  eval1.template gradients<0, false, false>(gradients_quad,
1652  values_dofs);
1653  if (integrate_val == true)
1654  eval1.template values<0, false, true>(values_quad,
1655  values_dofs);
1656  break;
1657  case 1:
1658  values_dofs[0] = values_quad[0];
1659  values_dofs[1] = gradients_quad[0];
1660  break;
1661  default:
1662  AssertThrow(false, ExcNotImplemented());
1663  }
1664  values_dofs += 2 * size_deg;
1665  values_quad += n_q_points;
1666  gradients_quad += dim * n_q_points;
1667  }
1668  }
1669  };
1670 
1671 
1672 
1673  template <int dim, int fe_degree, int n_components, typename Number>
1674  struct FEFaceNormalEvaluationImpl
1675  {
1676  template <bool do_evaluate, bool add_into_output>
1677  static void
1678  interpolate(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1679  const Number * input,
1680  Number * output,
1681  const bool do_gradients,
1682  const unsigned int face_no)
1683  {
1684  internal::EvaluatorTensorProduct<internal::evaluate_general,
1685  dim,
1686  fe_degree + 1,
1687  0,
1688  Number>
1689  evalf(data.shape_data_on_face[face_no % 2],
1690  AlignedVector<Number>(),
1691  AlignedVector<Number>(),
1692  data.fe_degree + 1,
1693  0);
1694 
1695  const unsigned int in_stride = do_evaluate ?
1696  data.dofs_per_component_on_cell :
1697  2 * data.dofs_per_component_on_face;
1698  const unsigned int out_stride = do_evaluate ?
1699  2 * data.dofs_per_component_on_face :
1700  data.dofs_per_component_on_cell;
1701  const unsigned int face_direction = face_no / 2;
1702  for (unsigned int c = 0; c < n_components; c++)
1703  {
1704  if (do_gradients)
1705  {
1706  if (face_direction == 0)
1707  evalf.template apply_face<0, do_evaluate, add_into_output, 1>(
1708  input, output);
1709  else if (face_direction == 1)
1710  evalf.template apply_face<1, do_evaluate, add_into_output, 1>(
1711  input, output);
1712  else
1713  evalf.template apply_face<2, do_evaluate, add_into_output, 1>(
1714  input, output);
1715  }
1716  else
1717  {
1718  if (face_direction == 0)
1719  evalf.template apply_face<0, do_evaluate, add_into_output, 0>(
1720  input, output);
1721  else if (face_direction == 1)
1722  evalf.template apply_face<1, do_evaluate, add_into_output, 0>(
1723  input, output);
1724  else
1725  evalf.template apply_face<2, do_evaluate, add_into_output, 0>(
1726  input, output);
1727  }
1728  input += in_stride;
1729  output += out_stride;
1730  }
1731  }
1732  };
1733 } // end of namespace internal
1734 
1735 
1736 DEAL_II_NAMESPACE_CLOSE
1737 
1738 #endif
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
internal::MatrixFreeFunctions::ElementType
ElementType
Definition: shape_info.h:50
internal::MatrixFreeFunctions::ShapeInfo::shape_values_eo
AlignedVector< Number > shape_values_eo
Definition: shape_info.h:180
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians
AlignedVector< Number > shape_hessians
Definition: shape_info.h:173
internal::FEEvaluationImplBasisChange::do_forward
static void do_forward(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:686
internal::MatrixFreeFunctions::truncated_tensor
Definition: shape_info.h:85
internal::MatrixFreeFunctions::ShapeInfo::shape_values
AlignedVector< Number > shape_values
Definition: shape_info.h:155
internal::FEEvaluationImplBasisChange::do_mass
static void do_mass(const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
Definition: evaluation_kernels.h:911
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1031
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::FEEvaluationImplBasisChange::do_backward
static void do_backward(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:805
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients_collocation_eo
AlignedVector< Number > shape_gradients_collocation_eo
Definition: shape_info.h:201
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
internal::MatrixFreeFunctions::tensor_symmetric
Definition: shape_info.h:74
Utilities::pow
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:428
internal::MatrixFreeFunctions::tensor_general
Definition: shape_info.h:79
Utilities
Definition: cuda.h:31
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians_eo
AlignedVector< Number > shape_hessians_eo
Definition: shape_info.h:194
internal::MatrixFreeFunctions::tensor_symmetric_collocation
Definition: shape_info.h:60
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients_eo
AlignedVector< Number > shape_gradients_eo
Definition: shape_info.h:187
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients
AlignedVector< Number > shape_gradients
Definition: shape_info.h:164
AlignedVector::size
size_type size() const
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians_collocation_eo
AlignedVector< Number > shape_hessians_collocation_eo
Definition: shape_info.h:208
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0
Definition: shape_info.h:92
Generated by   doxygen 1.8.14