Reference documentation for deal.II version 9.0.0
evaluation_kernels.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018 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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/vectorization.h>
22 #include <deal.II/base/utilities.h>
23 #include <deal.II/matrix_free/tensor_product_kernels.h>
24 #include <deal.II/matrix_free/shape_info.h>
25 
26 
27 DEAL_II_NAMESPACE_OPEN
28 
29 
30 
31 namespace internal
32 {
33  // Select evaluator type from element shape function type
34  template <MatrixFreeFunctions::ElementType element, bool is_long>
35  struct EvaluatorSelector {};
36 
37  template <bool is_long>
38  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general,is_long>
39  {
40  static const EvaluatorVariant variant = evaluate_general;
41  };
42 
43  template <>
44  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,false>
45  {
46  static const EvaluatorVariant variant = evaluate_symmetric;
47  };
48 
49  template <> struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,true>
50  {
51  static const EvaluatorVariant variant = evaluate_evenodd;
52  };
53 
54  template <bool is_long>
55  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor,is_long>
56  {
57  static const EvaluatorVariant variant = evaluate_general;
58  };
59 
60  template <>
61  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,false>
62  {
63  static const EvaluatorVariant variant = evaluate_general;
64  };
65 
66  template <>
67  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,true>
68  {
69  static const EvaluatorVariant variant = evaluate_evenodd;
70  };
71 
72  template <bool is_long>
73  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,is_long>
74  {
75  static const EvaluatorVariant variant = evaluate_evenodd;
76  };
77 
78 
79 
98  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
99  int n_q_points_1d, int n_components, typename Number>
100  struct FEEvaluationImpl
101  {
102  static
103  void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
104  const Number *values_dofs_actual,
105  Number *values_quad,
106  Number *gradients_quad,
107  Number *hessians_quad,
108  Number *scratch_data,
109  const bool evaluate_values,
110  const bool evaluate_gradients,
111  const bool evaluate_hessians);
112 
113  static
114  void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
115  Number *values_dofs_actual,
116  Number *values_quad,
117  Number *gradients_quad,
118  Number *scratch_data,
119  const bool integrate_values,
120  const bool integrate_gradients,
121  const bool add_into_values_array);
122  };
123 
124 
125 
126  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
127  int n_q_points_1d, int n_components, typename Number>
128  inline
129  void
130  FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
131  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
132  const Number *values_dofs_actual,
133  Number *values_quad,
134  Number *gradients_quad,
135  Number *hessians_quad,
136  Number *scratch_data,
137  const bool evaluate_values,
138  const bool evaluate_gradients,
139  const bool evaluate_hessians)
140  {
141  if (evaluate_values == false && evaluate_gradients == false && evaluate_hessians == false)
142  return;
143 
144  const EvaluatorVariant variant =
145  EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
146  typedef EvaluatorTensorProduct<variant, dim, fe_degree+1, n_q_points_1d,
147  Number> Eval;
148  Eval eval (variant == evaluate_evenodd ? shape_info.shape_values_eo :
149  shape_info.shape_values,
150  variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
151  shape_info.shape_gradients,
152  variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
153  shape_info.shape_hessians,
154  shape_info.fe_degree+1,
155  shape_info.n_q_points_1d);
156 
157  const unsigned int temp_size = Eval::n_rows_of_product == numbers::invalid_unsigned_int ? 0
158  : (Eval::n_rows_of_product > Eval::n_columns_of_product ?
159  Eval::n_rows_of_product : Eval::n_columns_of_product);
160  Number *temp1;
161  Number *temp2;
162  if (temp_size == 0)
163  {
164  temp1 = scratch_data;
165  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
166  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
167  }
168  else
169  {
170  temp1 = scratch_data;
171  temp2 = temp1 + temp_size;
172  }
173 
174  const unsigned int n_q_points = temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
175  const unsigned int dofs_per_comp = (type == MatrixFreeFunctions::truncated_tensor) ?
176  Utilities::fixed_power<dim>(shape_info.fe_degree+1) : shape_info.dofs_per_component_on_cell;
177  const Number *values_dofs = values_dofs_actual;
178  if (type == MatrixFreeFunctions::truncated_tensor)
179  {
180  Number *values_dofs_tmp = scratch_data+2*(std::max(shape_info.dofs_per_component_on_cell, shape_info.n_q_points));
181  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
182  unsigned int count_p = 0, count_q = 0;
183  for (int i=0; i<(dim>2?degree+1:1); ++i)
184  {
185  for (int j=0; j<(dim>1?degree+1-i:1); ++j)
186  {
187  for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
188  for (unsigned int c=0; c<n_components; ++c)
189  values_dofs_tmp[c*dofs_per_comp+count_q] = values_dofs_actual[c*shape_info.dofs_per_component_on_cell+count_p];
190  for (int k=degree+1-j-i; k<degree+1; ++k, ++count_q)
191  for (unsigned int c=0; c<n_components; ++c)
192  values_dofs_tmp[c*dofs_per_comp+count_q] = Number();
193  }
194  for (int j=degree+1-i; j<degree+1; ++j)
195  for (int k=0; k<degree+1; ++k, ++count_q)
196  for (unsigned int c=0; c<n_components; ++c)
197  values_dofs_tmp[c*dofs_per_comp+count_q] = Number();
198  }
199  AssertDimension(count_q, dofs_per_comp);
200  values_dofs = values_dofs_tmp;
201  }
202 
203  switch (dim)
204  {
205  case 1:
206  for (unsigned int c=0; c<n_components; c++)
207  {
208  if (evaluate_values == true)
209  eval.template values<0,true,false> (values_dofs, values_quad);
210  if (evaluate_gradients == true)
211  eval.template gradients<0,true,false>(values_dofs, gradients_quad);
212  if (evaluate_hessians == true)
213  eval.template hessians<0,true,false> (values_dofs, hessians_quad);
214 
215  // advance the next component in 1D array
216  values_dofs += dofs_per_comp;
217  values_quad += n_q_points;
218  gradients_quad += n_q_points;
219  hessians_quad += n_q_points;
220  }
221  break;
222 
223  case 2:
224  for (unsigned int c=0; c<n_components; c++)
225  {
226  // grad x
227  if (evaluate_gradients == true)
228  {
229  eval.template gradients<0,true,false> (values_dofs, temp1);
230  eval.template values<1,true,false> (temp1, gradients_quad);
231  }
232  if (evaluate_hessians == true)
233  {
234  // grad xy
235  if (evaluate_gradients == false)
236  eval.template gradients<0,true,false>(values_dofs, temp1);
237  eval.template gradients<1,true,false> (temp1, hessians_quad+2*n_q_points);
238 
239  // grad xx
240  eval.template hessians<0,true,false>(values_dofs, temp1);
241  eval.template values<1,true,false> (temp1, hessians_quad);
242  }
243 
244  // grad y
245  eval.template values<0,true,false> (values_dofs, temp1);
246  if (evaluate_gradients == true)
247  eval.template gradients<1,true,false> (temp1, gradients_quad+n_q_points);
248 
249  // grad yy
250  if (evaluate_hessians == true)
251  eval.template hessians<1,true,false> (temp1, hessians_quad+n_q_points);
252 
253  // val: can use values applied in x
254  if (evaluate_values == true)
255  eval.template values<1,true,false> (temp1, values_quad);
256 
257  // advance to the next component in 1D array
258  values_dofs += dofs_per_comp;
259  values_quad += n_q_points;
260  gradients_quad += 2*n_q_points;
261  hessians_quad += 3*n_q_points;
262  }
263  break;
264 
265  case 3:
266  for (unsigned int c=0; c<n_components; c++)
267  {
268  if (evaluate_gradients == true)
269  {
270  // grad x
271  eval.template gradients<0,true,false> (values_dofs, temp1);
272  eval.template values<1,true,false> (temp1, temp2);
273  eval.template values<2,true,false> (temp2, gradients_quad);
274  }
275 
276  if (evaluate_hessians == true)
277  {
278  // grad xz
279  if (evaluate_gradients == false)
280  {
281  eval.template gradients<0,true,false> (values_dofs, temp1);
282  eval.template values<1,true,false> (temp1, temp2);
283  }
284  eval.template gradients<2,true,false> (temp2, hessians_quad+4*n_q_points);
285 
286  // grad xy
287  eval.template gradients<1,true,false> (temp1, temp2);
288  eval.template values<2,true,false> (temp2, hessians_quad+3*n_q_points);
289 
290  // grad xx
291  eval.template hessians<0,true,false>(values_dofs, temp1);
292  eval.template values<1,true,false> (temp1, temp2);
293  eval.template values<2,true,false> (temp2, hessians_quad);
294  }
295 
296  // grad y
297  eval.template values<0,true,false> (values_dofs, temp1);
298  if (evaluate_gradients == true)
299  {
300  eval.template gradients<1,true,false>(temp1, temp2);
301  eval.template values<2,true,false> (temp2, gradients_quad+n_q_points);
302  }
303 
304  if (evaluate_hessians == true)
305  {
306  // grad yz
307  if (evaluate_gradients == false)
308  eval.template gradients<1,true,false>(temp1, temp2);
309  eval.template gradients<2,true,false> (temp2, hessians_quad+5*n_q_points);
310 
311  // grad yy
312  eval.template hessians<1,true,false> (temp1, temp2);
313  eval.template values<2,true,false> (temp2, hessians_quad+n_q_points);
314  }
315 
316  // grad z: can use the values applied in x direction stored in temp1
317  eval.template values<1,true,false> (temp1, temp2);
318  if (evaluate_gradients == true)
319  eval.template gradients<2,true,false> (temp2, gradients_quad+2*n_q_points);
320 
321  // grad zz: can use the values applied in x and y direction stored
322  // in temp2
323  if (evaluate_hessians == true)
324  eval.template hessians<2,true,false>(temp2, hessians_quad+2*n_q_points);
325 
326  // val: can use the values applied in x & y direction stored in temp2
327  if (evaluate_values == true)
328  eval.template values<2,true,false> (temp2, values_quad);
329 
330  // advance to the next component in 1D array
331  values_dofs += dofs_per_comp;
332  values_quad += n_q_points;
333  gradients_quad += 3*n_q_points;
334  hessians_quad += 6*n_q_points;
335  }
336  break;
337 
338  default:
339  AssertThrow(false, ExcNotImplemented());
340  }
341 
342  // case additional dof for FE_Q_DG0: add values; gradients and second
343  // derivatives evaluate to zero
344  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_values)
345  {
346  values_quad -= n_components*n_q_points;
347  values_dofs -= n_components*dofs_per_comp;
348  for (unsigned int c=0; c<n_components; ++c)
349  for (unsigned int q=0; q<shape_info.n_q_points; ++q)
350  values_quad[c*shape_info.n_q_points+q] +=
351  values_dofs[(c+1)*shape_info.dofs_per_component_on_cell-1];
352  }
353  }
354 
355 
356 
357  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
358  int n_q_points_1d, int n_components, typename Number>
359  inline
360  void
361  FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
362  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
363  Number *values_dofs_actual,
364  Number *values_quad,
365  Number *gradients_quad,
366  Number *scratch_data,
367  const bool integrate_values,
368  const bool integrate_gradients,
369  const bool add_into_values_array)
370  {
371  const EvaluatorVariant variant =
372  EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
373  typedef EvaluatorTensorProduct<variant, dim, fe_degree+1, n_q_points_1d,
374  Number> Eval;
375  Eval eval (variant == evaluate_evenodd ? shape_info.shape_values_eo :
376  shape_info.shape_values,
377  variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
378  shape_info.shape_gradients,
379  variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
380  shape_info.shape_hessians,
381  shape_info.fe_degree+1,
382  shape_info.n_q_points_1d);
383 
384  const unsigned int temp_size = Eval::n_rows_of_product == numbers::invalid_unsigned_int ? 0
385  : (Eval::n_rows_of_product > Eval::n_columns_of_product ?
386  Eval::n_rows_of_product : Eval::n_columns_of_product);
387  Number *temp1;
388  Number *temp2;
389  if (temp_size == 0)
390  {
391  temp1 = scratch_data;
392  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
393  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
394  }
395  else
396  {
397  temp1 = scratch_data;
398  temp2 = temp1 + temp_size;
399  }
400 
401  const unsigned int n_q_points = temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
402  const unsigned int dofs_per_comp = (type == MatrixFreeFunctions::truncated_tensor) ?
403  Utilities::fixed_power<dim>(shape_info.fe_degree+1) : shape_info.dofs_per_component_on_cell;
404  // expand dof_values to tensor product for truncated tensor products
405  Number *values_dofs = (type == MatrixFreeFunctions::truncated_tensor) ?
406  scratch_data+2*(std::max(shape_info.dofs_per_component_on_cell,
407  shape_info.n_q_points)) :
408  values_dofs_actual;
409 
410  switch (dim)
411  {
412  case 1:
413  for (unsigned int c=0; c<n_components; c++)
414  {
415  if (integrate_values == true)
416  {
417  if (add_into_values_array == false)
418  eval.template values<0,false,false> (values_quad, values_dofs);
419  else
420  eval.template values<0,false,true> (values_quad, values_dofs);
421  }
422  if (integrate_gradients == true)
423  {
424  if (integrate_values == true || add_into_values_array == true)
425  eval.template gradients<0,false,true> (gradients_quad, values_dofs);
426  else
427  eval.template gradients<0,false,false> (gradients_quad, values_dofs);
428  }
429 
430  // advance to the next component in 1D array
431  values_dofs += dofs_per_comp;
432  values_quad += n_q_points;
433  gradients_quad += n_q_points;
434  }
435  break;
436 
437  case 2:
438  for (unsigned int c=0; c<n_components; c++)
439  {
440  if (integrate_values == true &&
441  integrate_gradients == false)
442  {
443  eval.template values<1,false,false> (values_quad, temp1);
444  if (add_into_values_array == false)
445  eval.template values<0,false,false>(temp1, values_dofs);
446  else
447  eval.template values<0,false,true>(temp1, values_dofs);
448  }
449  if (integrate_gradients == true)
450  {
451  eval.template gradients<1,false,false> (gradients_quad+n_q_points, temp1);
452  if (integrate_values)
453  eval.template values<1,false,true> (values_quad, temp1);
454  if (add_into_values_array == false)
455  eval.template values<0,false,false>(temp1, values_dofs);
456  else
457  eval.template values<0,false,true>(temp1, values_dofs);
458  eval.template values<1,false,false> (gradients_quad, temp1);
459  eval.template gradients<0,false,true> (temp1, values_dofs);
460  }
461 
462  // advance to the next component in 1D array
463  values_dofs += dofs_per_comp;
464  values_quad += n_q_points;
465  gradients_quad += 2*n_q_points;
466  }
467  break;
468 
469  case 3:
470  for (unsigned int c=0; c<n_components; c++)
471  {
472  if (integrate_values == true &&
473  integrate_gradients == false)
474  {
475  eval.template values<2,false,false> (values_quad, temp1);
476  eval.template values<1,false,false> (temp1, temp2);
477  if (add_into_values_array == false)
478  eval.template values<0,false,false>(temp2, values_dofs);
479  else
480  eval.template values<0,false,true> (temp2, values_dofs);
481  }
482  if (integrate_gradients == true)
483  {
484  eval.template gradients<2,false,false>(gradients_quad+2*n_q_points, temp1);
485  if (integrate_values)
486  eval.template values<2,false,true> (values_quad, temp1);
487  eval.template values<1,false,false> (temp1, temp2);
488  eval.template values<2,false,false> (gradients_quad+n_q_points, temp1);
489  eval.template gradients<1,false,true> (temp1, temp2);
490  if (add_into_values_array == false)
491  eval.template values<0,false,false> (temp2, values_dofs);
492  else
493  eval.template values<0,false,true> (temp2, values_dofs);
494  eval.template values<2,false,false> (gradients_quad, temp1);
495  eval.template values<1,false,false> (temp1, temp2);
496  eval.template gradients<0,false,true> (temp2, values_dofs);
497  }
498 
499  // advance to the next component in 1D array
500  values_dofs += dofs_per_comp;
501  values_quad += n_q_points;
502  gradients_quad += 3*n_q_points;
503  }
504  break;
505 
506  default:
507  AssertThrow(false, ExcNotImplemented());
508  }
509 
510  // case FE_Q_DG0: add values, gradients and second derivatives are zero
511  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
512  {
513  values_dofs -= n_components * dofs_per_comp - shape_info.dofs_per_component_on_cell + 1;
514  values_quad -= n_components * n_q_points;
515  if (integrate_values)
516  for (unsigned int c=0; c<n_components; ++c)
517  {
518  values_dofs[0] = values_quad[0];
519  for (unsigned int q=1; q<shape_info.n_q_points; ++q)
520  values_dofs[0] += values_quad[q];
521  values_dofs += dofs_per_comp;
522  values_quad += n_q_points;
523  }
524  else
525  {
526  for (unsigned int c=0; c<n_components; ++c)
527  values_dofs[c*shape_info.dofs_per_component_on_cell] = Number();
528  values_dofs += n_components*shape_info.dofs_per_component_on_cell;
529  }
530  }
531 
532  if (type == MatrixFreeFunctions::truncated_tensor)
533  {
534  values_dofs -= dofs_per_comp*n_components;
535  unsigned int count_p = 0, count_q = 0;
536  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
537  for (int i=0; i<(dim>2?degree+1:1); ++i)
538  {
539  for (int j=0; j<(dim>1?degree+1-i:1); ++j)
540  {
541  for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
542  {
543  for (unsigned int c=0; c<n_components; ++c)
544  values_dofs_actual[c*shape_info.dofs_per_component_on_cell+count_p] = values_dofs[c*dofs_per_comp+count_q];
545  }
546  count_q += j+i;
547  }
548  count_q += i*(degree+1);
549  }
550  AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
551  }
552  }
553 
554 
555 
567  template <EvaluatorVariant variant, int dim, int basis_size_1, int basis_size_2, int n_components,
568  typename Number, typename Number2>
569  struct FEEvaluationImplBasisChange
570  {
571  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
572  "The second dimension must not be smaller than the first");
573 
596 #ifndef DEBUG
597  DEAL_II_ALWAYS_INLINE
598 #endif
599  static void do_forward (const AlignedVector<Number2> &transformation_matrix,
600  const Number *values_in,
601  Number *values_out,
602  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
603  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
604  {
605  Assert(basis_size_1 != 0 ||
606  basis_size_1_variable <= basis_size_2_variable,
607  ExcMessage("The second dimension must not be smaller than the first"));
608 
609  // we do recursion until dim==1 or dim==2 and we have
610  // basis_size_1==basis_size_2. The latter optimization increases
611  // optimization possibilities for the compiler but does only work for
612  // aliased pointers if the sizes are equal.
613  constexpr int next_dim = (dim > 2 || ((basis_size_1 == 0 || basis_size_2>basis_size_1)
614  && dim>1)) ? dim-1 : dim;
615 
616  EvaluatorTensorProduct<variant, dim, basis_size_1, (basis_size_1==0 ? 0 : basis_size_2),
617  Number,Number2> eval_val (transformation_matrix,
618  AlignedVector<Number2>(),
619  AlignedVector<Number2>(),
620  basis_size_1_variable,
621  basis_size_2_variable);
622  const unsigned int np_1 = basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
623  const unsigned int np_2 = basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
624  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
625  ExcMessage("Cannot transform with 0-point basis"));
626  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
627  ExcMessage("Cannot transform with 0-point basis"));
628 
629  // run loop backwards to ensure correctness if values_in aliases with
630  // values_out in case with basis_size_1 < basis_size_2
631  values_in = values_in + n_components*Utilities::fixed_power<dim>(np_1);
632  values_out = values_out + n_components*Utilities::fixed_power<dim>(np_2);
633  for (unsigned int c=n_components; c!=0; --c)
634  {
635  values_in -= Utilities::fixed_power<dim>(np_1);
636  values_out -= Utilities::fixed_power<dim>(np_2);
637  if (next_dim < dim)
638  for (unsigned int q=np_1; q!=0; --q)
639  FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,1,Number,Number2>
640  ::do_forward(transformation_matrix,
641  values_in + (q-1)*Utilities::fixed_power<next_dim>(np_1),
642  values_out + (q-1)*Utilities::fixed_power<next_dim>(np_2),
643  basis_size_1_variable,
644  basis_size_2_variable);
645 
646  // the recursion stops if dim==1 or if dim==2 and
647  // basis_size_1==basis_size_2 (the latter is used because the
648  // compiler generates nicer code)
649  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
650  {
651  eval_val.template values<0,true,false>(values_in, values_out);
652  eval_val.template values<1,true,false>(values_out, values_out);
653  }
654  else if (dim==1)
655  eval_val.template values<dim-1,true,false>(values_in, values_out);
656  else
657  eval_val.template values<dim-1,true,false>(values_out, values_out);
658  }
659  }
660 
691 #ifndef DEBUG
692  DEAL_II_ALWAYS_INLINE
693 #endif
694  static void do_backward (const AlignedVector<Number2> &transformation_matrix,
695  const bool add_into_result,
696  Number *values_in,
697  Number *values_out,
698  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
699  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
700  {
701  Assert(basis_size_1 != 0 ||
702  basis_size_1_variable <= basis_size_2_variable,
703  ExcMessage("The second dimension must not be smaller than the first"));
704  Assert(add_into_result == false || values_in != values_out,
705  ExcMessage("Input and output cannot alias with each other when "
706  "adding the result of the basis change to existing data"));
707 
708  constexpr int next_dim = (dim > 2 || ((basis_size_1 == 0 || basis_size_2>basis_size_1)
709  && dim>1)) ? dim-1 : dim;
710  EvaluatorTensorProduct<variant, dim, basis_size_1, (basis_size_1==0 ? 0 : basis_size_2),
711  Number,Number2> eval_val (transformation_matrix,
712  AlignedVector<Number2>(),
713  AlignedVector<Number2>(),
714  basis_size_1_variable,
715  basis_size_2_variable);
716  const unsigned int np_1 = basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
717  const unsigned int np_2 = basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
718  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
719  ExcMessage("Cannot transform with 0-point basis"));
720  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
721  ExcMessage("Cannot transform with 0-point basis"));
722 
723  for (unsigned int c=0; c<n_components; ++c)
724  {
725  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
726  {
727  eval_val.template values<1,false,false>(values_in, values_in);
728  if (add_into_result)
729  eval_val.template values<0,false,true>(values_in, values_out);
730  else
731  eval_val.template values<0,false,false>(values_in, values_out);
732  }
733  else
734  {
735  if (dim==1 && add_into_result)
736  eval_val.template values<0,false,true>(values_in, values_out);
737  else if (dim==1)
738  eval_val.template values<0,false,false>(values_in, values_out);
739  else
740  eval_val.template values<dim-1,false,false>(values_in, values_in);
741  }
742  if (next_dim < dim)
743  for (unsigned int q=0; q<np_1; ++q)
744  FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,1,Number,Number2>
745  ::do_backward(transformation_matrix,
746  add_into_result,
747  values_in + q*Utilities::fixed_power<next_dim>(np_2),
748  values_out + q*Utilities::fixed_power<next_dim>(np_1),
749  basis_size_1_variable, basis_size_2_variable);
750 
751  values_in += Utilities::fixed_power<dim>(np_2);
752  values_out += Utilities::fixed_power<dim>(np_1);
753  }
754  }
755 
775  static void do_mass (const AlignedVector<Number2> &transformation_matrix,
776  const AlignedVector<Number> &coefficients,
777  const Number *values_in,
778  Number *scratch_data,
779  Number *values_out)
780  {
781  constexpr int next_dim = dim > 1 ? dim-1 : dim;
782  Number *my_scratch = basis_size_1 != basis_size_2 ? scratch_data : values_out;
783  for (unsigned int q=basis_size_1; q!=0; --q)
784  FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,n_components,Number,Number2>
785  ::do_forward(transformation_matrix,
786  values_in + (q-1)*Utilities::pow(basis_size_1, dim-1),
787  my_scratch + (q-1)*Utilities::pow(basis_size_2, dim-1));
788  EvaluatorTensorProduct<variant, dim, basis_size_1, basis_size_2,
789  Number,Number2> eval_val (transformation_matrix);
790  const unsigned int n_inner_blocks = (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
791  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim-1);
792  for (unsigned int ii=0; ii<n_blocks; ii+=n_inner_blocks)
793  for (unsigned int c=0; c<n_components; ++c)
794  {
795  for (unsigned int i=ii; i<ii+n_inner_blocks; ++i)
796  eval_val.template values_one_line<dim-1,true,false> (my_scratch+i, my_scratch+i);
797  for (unsigned int q=0; q<basis_size_2; ++q)
798  for (unsigned int i=ii; i<ii+n_inner_blocks; ++i)
799  my_scratch[i+q*n_blocks] *= coefficients[i+q*n_blocks];
800  for (unsigned int i=ii; i<ii+n_inner_blocks; ++i)
801  eval_val.template values_one_line<dim-1,false,false>(my_scratch+i, my_scratch+i);
802  }
803  for (unsigned int q=0; q<basis_size_1; ++q)
804  FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,n_components,Number,Number2>
805  ::do_backward(transformation_matrix, false,
806  my_scratch + q*Utilities::pow(basis_size_2, dim-1),
807  values_out + q*Utilities::pow(basis_size_1, dim-1));
808  }
809  };
810 
811 
812 
827  template <int dim, int fe_degree, int n_components, typename Number>
828  struct FEEvaluationImplCollocation
829  {
830  static
831  void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
832  const Number *values_dofs,
833  Number *values_quad,
834  Number *gradients_quad,
835  Number *hessians_quad,
836  Number *scratch_data,
837  const bool evaluate_values,
838  const bool evaluate_gradients,
839  const bool evaluate_hessians);
840 
841  static
842  void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
843  Number *values_dofs,
844  Number *values_quad,
845  Number *gradients_quad,
846  Number *scratch_data,
847  const bool integrate_values,
848  const bool integrate_gradients,
849  const bool add_into_values_array);
850  };
851 
852 
853 
854  template <int dim, int fe_degree, int n_components, typename Number>
855  inline
856  void
857  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
858  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
859  const Number *values_dofs,
860  Number *values_quad,
861  Number *gradients_quad,
862  Number *hessians_quad,
863  Number *,
864  const bool evaluate_values,
865  const bool evaluate_gradients,
866  const bool evaluate_hessians)
867  {
868  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
869  (fe_degree+2)/2*(fe_degree+1));
870 
871  EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree+1, fe_degree+1, Number>
872  eval(AlignedVector<Number>(),
873  shape_info.shape_gradients_collocation_eo,
874  shape_info.shape_hessians_collocation_eo);
875  constexpr unsigned int n_q_points = Utilities::pow(fe_degree+1, dim);
876 
877  for (unsigned int c=0; c<n_components; c++)
878  {
879  if (evaluate_values == true)
880  for (unsigned int i=0; i<n_q_points; ++i)
881  values_quad[i] = values_dofs[i];
882  if (evaluate_gradients == true || evaluate_hessians == true)
883  {
884  eval.template gradients<0,true,false>(values_dofs, gradients_quad);
885  if (dim > 1)
886  eval.template gradients<1,true,false>(values_dofs, gradients_quad+n_q_points);
887  if (dim > 2)
888  eval.template gradients<2,true,false>(values_dofs, gradients_quad+2*n_q_points);
889  }
890  if (evaluate_hessians == true)
891  {
892  eval.template hessians<0,true,false> (values_dofs, hessians_quad);
893  if (dim > 1)
894  {
895  eval.template gradients<1,true,false> (gradients_quad, hessians_quad+dim*n_q_points);
896  eval.template hessians<1,true,false> (values_dofs, hessians_quad+n_q_points);
897  }
898  if (dim > 2)
899  {
900  eval.template gradients<2,true,false> (gradients_quad, hessians_quad+4*n_q_points);
901  eval.template gradients<2,true,false> (gradients_quad+n_q_points, hessians_quad+5*n_q_points);
902  eval.template hessians<2,true,false> (values_dofs, hessians_quad+2*n_q_points);
903  }
904  hessians_quad += (dim*(dim+1))/2*n_q_points;
905  }
906  gradients_quad += dim*n_q_points;
907  values_quad += n_q_points;
908  values_dofs += n_q_points;
909  }
910  }
911 
912 
913 
914  template <int dim, int fe_degree, int n_components, typename Number>
915  inline
916  void
917  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
918  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
919  Number *values_dofs,
920  Number *values_quad,
921  Number *gradients_quad,
922  Number *,
923  const bool integrate_values,
924  const bool integrate_gradients,
925  const bool add_into_values_array)
926  {
927  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
928  (fe_degree+2)/2*(fe_degree+1));
929 
930  EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree+1, fe_degree+1, Number>
931  eval(AlignedVector<Number>(),
932  shape_info.shape_gradients_collocation_eo,
933  shape_info.shape_hessians_collocation_eo);
934  constexpr unsigned int n_q_points = Utilities::pow(fe_degree+1, dim);
935 
936  for (unsigned int c=0; c<n_components; c++)
937  {
938  if (integrate_values == true && add_into_values_array == false)
939  for (unsigned int i=0; i<n_q_points; ++i)
940  values_dofs[i] = values_quad[i];
941  else if (integrate_values == true)
942  for (unsigned int i=0; i<n_q_points; ++i)
943  values_dofs[i] += values_quad[i];
944  if (integrate_gradients == true)
945  {
946  if (integrate_values == true || add_into_values_array == true)
947  eval.template gradients<0,false,true>(gradients_quad, values_dofs);
948  else
949  eval.template gradients<0,false,false>(gradients_quad, values_dofs);
950  if (dim > 1)
951  eval.template gradients<1,false,true>(gradients_quad+n_q_points, values_dofs);
952  if (dim > 2)
953  eval.template gradients<2,false,true>(gradients_quad+2*n_q_points, values_dofs);
954  }
955  gradients_quad += dim*n_q_points;
956  values_quad += n_q_points;
957  values_dofs += n_q_points;
958  }
959  }
960 
961 
962 
963 
976  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
977  struct FEEvaluationImplTransformToCollocation
978  {
979  static
980  void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
981  const Number *values_dofs,
982  Number *values_quad,
983  Number *gradients_quad,
984  Number *hessians_quad,
985  Number *scratch_data,
986  const bool evaluate_values,
987  const bool evaluate_gradients,
988  const bool evaluate_hessians);
989 
990  static
991  void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
992  Number *values_dofs,
993  Number *values_quad,
994  Number *gradients_quad,
995  Number *scratch_data,
996  const bool integrate_values,
997  const bool integrate_gradients,
998  const bool add_into_values_array);
999  };
1000 
1001 
1002 
1003  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
1004  inline
1005  void
1006  FEEvaluationImplTransformToCollocation<dim, fe_degree, n_q_points_1d, n_components, Number>
1007  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1008  const Number *values_dofs,
1009  Number *values_quad,
1010  Number *gradients_quad,
1011  Number *hessians_quad,
1012  Number *,
1013  const bool ,
1014  const bool evaluate_gradients,
1015  const bool evaluate_hessians)
1016  {
1017  Assert(n_q_points_1d > fe_degree,
1018  ExcMessage("You lose information when going to a collocation space "
1019  "of lower degree, so the evaluation results would be "
1020  "wrong. Thus, this class does not permit the desired "
1021  "operation."));
1022  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1023 
1024  for (unsigned int c=0; c<n_components; c++)
1025  {
1026  FEEvaluationImplBasisChange<evaluate_evenodd, dim,
1027  (fe_degree>=n_q_points_1d?n_q_points_1d:fe_degree+1),
1028  n_q_points_1d,1,Number,Number>
1029  ::do_forward(shape_info.shape_values_eo,
1030  values_dofs, values_quad);
1031 
1032  // apply derivatives in the collocation space
1033  if (evaluate_gradients == true || evaluate_hessians == true)
1034  FEEvaluationImplCollocation<dim,n_q_points_1d-1,1,Number>::
1035  evaluate(shape_info, values_quad, nullptr, gradients_quad, hessians_quad,
1036  nullptr, false, evaluate_gradients, evaluate_hessians);
1037 
1038  values_dofs += shape_info.dofs_per_component_on_cell;
1039  values_quad += n_q_points;
1040  gradients_quad += dim*n_q_points;
1041  hessians_quad += (dim*(dim+1))/2*n_q_points;
1042  }
1043  }
1044 
1045 
1046 
1047  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
1048  inline
1049  void
1050  FEEvaluationImplTransformToCollocation<dim, fe_degree, n_q_points_1d, n_components, Number>
1051  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1052  Number *values_dofs,
1053  Number *values_quad,
1054  Number *gradients_quad,
1055  Number *,
1056  const bool integrate_values,
1057  const bool integrate_gradients,
1058  const bool add_into_values_array)
1059  {
1060  Assert(n_q_points_1d > fe_degree,
1061  ExcMessage("You lose information when going to a collocation space "
1062  "of lower degree, so the evaluation results would be "
1063  "wrong. Thus, this class does not permit the desired "
1064  "operation."));
1065  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1066  (n_q_points_1d+1)/2*n_q_points_1d);
1067  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1068 
1069  for (unsigned int c=0; c<n_components; c++)
1070  {
1071 
1072  // apply derivatives in collocation space
1073  if (integrate_gradients == true)
1074  FEEvaluationImplCollocation<dim,n_q_points_1d-1,1,Number>::
1075  integrate(shape_info, values_quad, nullptr, gradients_quad, nullptr, false,
1076  integrate_gradients,/*add_into_values_array=*/integrate_values);
1077 
1078  // transform back to the original space
1079  FEEvaluationImplBasisChange<evaluate_evenodd, dim,
1080  (fe_degree>=n_q_points_1d?n_q_points_1d:fe_degree+1),
1081  n_q_points_1d,1,Number,Number>
1082  ::do_backward(shape_info.shape_values_eo,
1083  add_into_values_array,
1084  values_quad,
1085  values_dofs);
1086  gradients_quad += dim*n_q_points;
1087  values_quad += n_q_points;
1088  values_dofs += shape_info.dofs_per_component_on_cell;
1089  }
1090  }
1091 
1092 
1093 
1094  template <bool symmetric_evaluate, int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
1095  struct FEFaceEvaluationImpl
1096  {
1097  static
1098  void evaluate_in_face (const MatrixFreeFunctions::ShapeInfo<Number> &data,
1099  Number *values_dofs,
1100  Number *values_quad,
1101  Number *gradients_quad,
1102  Number *scratch_data,
1103  const bool evaluate_val,
1104  const bool evaluate_grad,
1105  const unsigned int subface_index)
1106  {
1107  const AlignedVector<Number> &val1
1108  = symmetric_evaluate ? data.shape_values_eo :
1109  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1110  data.shape_values : data.values_within_subface[subface_index%2]);
1111  const AlignedVector<Number> &val2
1112  = symmetric_evaluate ? data.shape_values_eo :
1113  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1114  data.shape_values : data.values_within_subface[subface_index/2]);
1115 
1116  const AlignedVector<Number> &grad1
1117  = symmetric_evaluate ? data.shape_gradients_eo :
1118  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1119  data.shape_gradients : data.gradients_within_subface[subface_index%2]);
1120  const AlignedVector<Number> &grad2
1121  = symmetric_evaluate ? data.shape_gradients_eo :
1122  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1123  data.shape_gradients : data.gradients_within_subface[subface_index/2]);
1124 
1125  typedef internal::EvaluatorTensorProduct
1126  <symmetric_evaluate ? internal::evaluate_evenodd :internal::evaluate_general,
1127  dim-1,fe_degree+1,n_q_points_1d,Number> Eval;
1128  Eval eval1(val1,grad1,AlignedVector<Number>(),
1129  data.fe_degree+1, data.n_q_points_1d);
1130  Eval eval2(val2,grad2,AlignedVector<Number>(),
1131  data.fe_degree+1, data.n_q_points_1d);
1132 
1133  const unsigned int size_deg = fe_degree > -1 ?
1134  Utilities::pow(fe_degree+1, dim-1) :
1135  (dim > 1 ? Utilities::fixed_power<dim-1>(data.fe_degree+1) : 1);
1136 
1137  const unsigned int n_q_points = fe_degree > -1 ?
1138  Utilities::pow(n_q_points_1d, dim-1) : data.n_q_points_face;
1139 
1140  if (evaluate_grad == false)
1141  for (unsigned int c=0; c<n_components; ++c)
1142  {
1143  switch (dim)
1144  {
1145  case 3:
1146  eval1.template values<0,true,false>(values_dofs, values_quad);
1147  eval2.template values<1,true,false>(values_quad, values_quad);
1148  break;
1149  case 2:
1150  eval1.template values<0,true,false>(values_dofs, values_quad);
1151  break;
1152  case 1:
1153  values_quad[c] = values_dofs[2*c];
1154  break;
1155  default:
1156  Assert(false, ExcNotImplemented());
1157  }
1158  values_dofs += 2*size_deg;
1159  values_quad += n_q_points;
1160  }
1161  else
1162  for (unsigned int c=0; c<n_components; ++c)
1163  {
1164  switch (dim)
1165  {
1166  case 3:
1167  if (symmetric_evaluate && n_q_points_1d > fe_degree)
1168  {
1169  eval1.template values<0,true,false>(values_dofs, values_quad);
1170  eval1.template values<1,true,false>(values_quad, values_quad);
1171  internal::EvaluatorTensorProduct
1172  <internal::evaluate_evenodd,dim-1,n_q_points_1d,n_q_points_1d,Number> eval_grad
1173  (AlignedVector<Number>(),
1174  data.shape_gradients_collocation_eo,
1175  AlignedVector<Number>());
1176  eval_grad.template gradients<0,true,false>(values_quad, gradients_quad);
1177  eval_grad.template gradients<1,true,false>(values_quad,
1178  gradients_quad+n_q_points);
1179  }
1180  else
1181  {
1182  eval1.template gradients<0,true,false>(values_dofs, scratch_data);
1183  eval2.template values<1,true,false>(scratch_data, gradients_quad);
1184 
1185  eval1.template values<0,true,false>(values_dofs, scratch_data);
1186  eval2.template gradients<1,true,false>(scratch_data, gradients_quad+n_q_points);
1187 
1188  if (evaluate_val == true)
1189  eval2.template values<1,true,false>(scratch_data, values_quad);
1190  }
1191  eval1.template values<0,true,false>(values_dofs+size_deg, scratch_data);
1192  eval2.template values<1,true,false>(scratch_data,
1193  gradients_quad+(dim-1)*n_q_points);
1194 
1195  break;
1196  case 2:
1197  eval1.template values<0,true,false>(values_dofs+size_deg,
1198  gradients_quad+(dim-1)*n_q_points);
1199  eval1.template gradients<0,true,false>(values_dofs, gradients_quad);
1200  if (evaluate_val == true)
1201  eval1.template values<0,true,false>(values_dofs, values_quad);
1202  break;
1203  case 1:
1204  values_quad[0] = values_dofs[0];
1205  gradients_quad[0] = values_dofs[1];
1206  break;
1207  default:
1208  AssertThrow(false, ExcNotImplemented());
1209  }
1210  values_dofs += 2*size_deg;
1211  values_quad += n_q_points;
1212  gradients_quad += dim*n_q_points;
1213  }
1214  }
1215 
1216  static
1217  void integrate_in_face (const MatrixFreeFunctions::ShapeInfo<Number> &data,
1218  Number *values_dofs,
1219  Number *values_quad,
1220  Number *gradients_quad,
1221  Number *scratch_data,
1222  const bool integrate_val,
1223  const bool integrate_grad,
1224  const unsigned int subface_index)
1225  {
1226  const AlignedVector<Number> &val1
1227  = symmetric_evaluate ? data.shape_values_eo :
1228  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1229  data.shape_values : data.values_within_subface[subface_index%2]);
1230  const AlignedVector<Number> &val2
1231  = symmetric_evaluate ? data.shape_values_eo :
1232  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1233  data.shape_values : data.values_within_subface[subface_index/2]);
1234 
1235  const AlignedVector<Number> &grad1
1236  = symmetric_evaluate ? data.shape_gradients_eo :
1237  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1238  data.shape_gradients : data.gradients_within_subface[subface_index%2]);
1239  const AlignedVector<Number> &grad2
1240  = symmetric_evaluate ? data.shape_gradients_eo :
1241  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1242  data.shape_gradients : data.gradients_within_subface[subface_index/2]);
1243 
1244  typedef internal::EvaluatorTensorProduct
1245  <symmetric_evaluate ? internal::evaluate_evenodd :internal::evaluate_general,
1246  dim-1,fe_degree+1,n_q_points_1d,Number> Eval;
1247  Eval eval1(val1,grad1,val1,data.fe_degree+1, data.n_q_points_1d);
1248  Eval eval2(val2,grad2,val1,data.fe_degree+1, data.n_q_points_1d);
1249 
1250  const unsigned int size_deg = fe_degree > -1 ?
1251  Utilities::pow(fe_degree+1, dim-1) :
1252  (dim > 1 ? Utilities::fixed_power<dim-1>(data.fe_degree+1) : 1);
1253 
1254  const unsigned int n_q_points = fe_degree > -1 ?
1255  Utilities::fixed_int_power<n_q_points_1d,dim-1>::value : data.n_q_points_face;
1256 
1257  if (integrate_grad == false)
1258  for (unsigned int c=0; c<n_components; ++c)
1259  {
1260  switch (dim)
1261  {
1262  case 3:
1263  eval2.template values<1,false,false>(values_quad, values_quad);
1264  eval1.template values<0,false,false>(values_quad, values_dofs);
1265  break;
1266  case 2:
1267  eval1.template values<0,false,false>(values_quad, values_dofs);
1268  break;
1269  case 1:
1270  values_dofs[2*c] = values_quad[c][0];
1271  break;
1272  default:
1273  Assert(false, ExcNotImplemented());
1274  }
1275  values_dofs += 2*size_deg;
1276  values_quad += n_q_points;
1277  }
1278  else
1279  for (unsigned int c=0; c<n_components; ++c)
1280  {
1281  switch (dim)
1282  {
1283  case 3:
1284  eval2.template values<1,false,false> (gradients_quad+2*n_q_points,
1285  gradients_quad+2*n_q_points);
1286  eval1.template values<0,false,false> (gradients_quad+2*n_q_points,
1287  values_dofs+size_deg);
1288  if (symmetric_evaluate && n_q_points_1d > fe_degree)
1289  {
1290  internal::EvaluatorTensorProduct <internal::evaluate_evenodd,
1291  dim-1,n_q_points_1d,n_q_points_1d,Number> eval_grad
1292  (AlignedVector<Number>(),
1293  data.shape_gradients_collocation_eo,
1294  AlignedVector<Number>());
1295  if (integrate_val)
1296  eval_grad.template gradients<1,false,true>(gradients_quad+n_q_points,
1297  values_quad);
1298  else
1299  eval_grad.template gradients<1,false,false>(gradients_quad+n_q_points,
1300  values_quad);
1301  eval_grad.template gradients<0,false,true>(gradients_quad,
1302  values_quad);
1303  eval1.template values<1,false,false>(values_quad, values_quad);
1304  eval1.template values<0,false,false>(values_quad, values_dofs);
1305  }
1306  else
1307  {
1308  if (integrate_val)
1309  {
1310  eval2.template values<1,false,false> (values_quad, scratch_data);
1311  eval2.template gradients<1,false,true> (gradients_quad+n_q_points,
1312  scratch_data);
1313  }
1314  else
1315  eval2.template gradients<1,false,false> (gradients_quad+n_q_points,
1316  scratch_data);
1317 
1318  eval1.template values<0,false,false> (scratch_data, values_dofs);
1319  eval2.template values<1,false,false> (gradients_quad, scratch_data);
1320  eval1.template gradients<0,false,true> (scratch_data, values_dofs);
1321  }
1322  break;
1323  case 2:
1324  eval1.template values<0,false,false>(gradients_quad+n_q_points,
1325  values_dofs+size_deg);
1326  eval1.template gradients<0,false,false>(gradients_quad, values_dofs);
1327  if (integrate_val == true)
1328  eval1.template values<0,false,true>(values_quad, values_dofs);
1329  break;
1330  case 1:
1331  values_dofs[0] = values_quad[0];
1332  values_dofs[1] = gradients_quad[0];
1333  break;
1334  default:
1335  AssertThrow(false, ExcNotImplemented());
1336  }
1337  values_dofs += 2*size_deg;
1338  values_quad += n_q_points;
1339  gradients_quad += dim*n_q_points;
1340  }
1341  }
1342  };
1343 
1344 
1345 
1346  template <int dim, int fe_degree, int n_components, typename Number>
1347  struct FEFaceNormalEvaluationImpl
1348  {
1349  template <bool do_evaluate, bool add_into_output>
1350  static void interpolate(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1351  const Number *input,
1352  Number *output,
1353  const bool do_gradients,
1354  const unsigned int face_no)
1355  {
1356  internal::EvaluatorTensorProduct<internal::evaluate_general,dim,
1357  fe_degree+1,0,Number>
1358  evalf(data.shape_data_on_face[face_no%2],
1359  AlignedVector<Number>(),
1360  AlignedVector<Number>(),
1361  data.fe_degree+1, 0);
1362 
1363  const unsigned int in_stride = do_evaluate ? data.dofs_per_component_on_cell : 2*data.dofs_per_component_on_face;
1364  const unsigned int out_stride = do_evaluate ? 2*data.dofs_per_component_on_face : data.dofs_per_component_on_cell;
1365  const unsigned int face_direction = face_no / 2;
1366  for (unsigned int c=0; c<n_components; c++)
1367  {
1368  if (do_gradients)
1369  {
1370  if (face_direction == 0)
1371  evalf.template apply_face<0,do_evaluate,add_into_output,1>(input, output);
1372  else if (face_direction == 1)
1373  evalf.template apply_face<1,do_evaluate,add_into_output,1>(input, output);
1374  else
1375  evalf.template apply_face<2,do_evaluate,add_into_output,1>(input, output);
1376  }
1377  else
1378  {
1379  if (face_direction == 0)
1380  evalf.template apply_face<0,do_evaluate,add_into_output,0>(input, output);
1381  else if (face_direction == 1)
1382  evalf.template apply_face<1,do_evaluate,add_into_output,0>(input, output);
1383  else
1384  evalf.template apply_face<2,do_evaluate,add_into_output,0>(input, output);
1385  }
1386  input += in_stride;
1387  output += out_stride;
1388  }
1389  }
1390  };
1391 } // end of namespace internal
1392 
1393 
1394 DEAL_II_NAMESPACE_CLOSE
1395 
1396 #endif
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
Utilities::fixed_int_power
Definition: utilities.h:330
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1248
Utilities::pow
constexpr unsigned int pow(const unsigned int base, const unsigned int iexp)
Definition: utilities.h:354
internal::MatrixFreeFunctions::ShapeInfo::shape_values_eo
AlignedVector< Number > shape_values_eo
Definition: shape_info.h:161
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1221
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians
AlignedVector< Number > shape_hessians
Definition: shape_info.h:154
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:599
internal::MatrixFreeFunctions::ShapeInfo::shape_values
AlignedVector< Number > shape_values
Definition: shape_info.h:136
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:775
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:863
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:694
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1142
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients_collocation_eo
AlignedVector< Number > shape_gradients_collocation_eo
Definition: shape_info.h:182
Utilities
Definition: cuda.h:31
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians_eo
AlignedVector< Number > shape_hessians_eo
Definition: shape_info.h:175
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients_eo
AlignedVector< Number > shape_gradients_eo
Definition: shape_info.h:168
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients
AlignedVector< Number > shape_gradients
Definition: shape_info.h:145
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:189
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)
Generated by   doxygen 1.8.14