Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 8.5.1
All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
evaluation_kernels.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 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_gausslobatto,is_long>
74  {
75  static const EvaluatorVariant variant = evaluate_evenodd;
76  };
77 
78 
79 
80  // This struct performs the evaluation of function values, gradients and
81  // Hessians for tensor-product finite elements. The operation is used for
82  // both the symmetric and non-symmetric case, which use different apply
83  // functions 'values', 'gradients' in the individual coordinate
84  // directions. The apply functions for values are provided through one of
85  // the template classes EvaluatorTensorProduct which in turn are selected
86  // from the MatrixFreeFunctions::ElementType template argument.
87  //
88  // There is a specialization made for Gauss-Lobatto elements further down
89  // where the 'values' operation is identity, which allows us to write
90  // shorter code.
91  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
92  int n_q_points_1d, int n_components, typename Number>
93  struct FEEvaluationImpl
94  {
95  static
96  void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
97  VectorizedArray<Number> *values_dofs_actual[],
98  VectorizedArray<Number> *values_quad[],
99  VectorizedArray<Number> *gradients_quad[][dim],
100  VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
101  VectorizedArray<Number> *scratch_data,
102  const bool evaluate_val,
103  const bool evaluate_grad,
104  const bool evaluate_lapl);
105 
106  static
107  void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
108  VectorizedArray<Number> *values_dofs_actual[],
109  VectorizedArray<Number> *values_quad[],
110  VectorizedArray<Number> *gradients_quad[][dim],
111  VectorizedArray<Number> *scratch_data,
112  const bool evaluate_val,
113  const bool evaluate_grad);
114  };
115 
116 
117  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
118  int n_q_points_1d, int n_components, typename Number>
119  inline
120  void
121  FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
122  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
123  VectorizedArray<Number> *values_dofs_actual[],
124  VectorizedArray<Number> *values_quad[],
125  VectorizedArray<Number> *gradients_quad[][dim],
126  VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
127  VectorizedArray<Number> *scratch_data,
128  const bool evaluate_val,
129  const bool evaluate_grad,
130  const bool evaluate_lapl)
131  {
132  if (evaluate_val == false && evaluate_grad == false && evaluate_lapl == false)
133  return;
134 
135  const EvaluatorVariant variant =
136  EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
137  typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
139  Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
140  shape_info.shape_values,
141  variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
142  shape_info.shape_gradients,
143  variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
144  shape_info.shape_hessians,
145  shape_info.fe_degree,
146  shape_info.n_q_points_1d);
147 
148  const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
149  : (Eval::dofs_per_cell > Eval::n_q_points ?
150  Eval::dofs_per_cell : Eval::n_q_points);
151  const unsigned int max_stack_size = 100;
152  VectorizedArray<Number> temp_data[(temp_size > 0 && temp_size < max_stack_size) ? 2*temp_size : 1];
155  if (temp_size == 0)
156  {
157  temp1 = scratch_data;
158  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
159  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
160  }
161  else if (temp_size < max_stack_size)
162  {
163  temp1 = &temp_data[0];
164  temp2 = temp1 + temp_size;
165  }
166  else
167  {
168  temp1 = scratch_data;
169  temp2 = temp1 + temp_size;
170  }
171 
172  VectorizedArray<Number> **values_dofs = values_dofs_actual;
173  VectorizedArray<Number> *expanded_dof_values[n_components];
174  if (type == MatrixFreeFunctions::truncated_tensor)
175  {
176  values_dofs = expanded_dof_values;
177  for (unsigned int c=0; c<n_components; ++c)
178  expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
179  shape_info.n_q_points)) +
180  c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
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  expanded_dof_values[c][count_q] = values_dofs_actual[c][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  expanded_dof_values[c][count_q] = VectorizedArray<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  expanded_dof_values[c][count_q] = VectorizedArray<Number>();
198  }
199  AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
200  }
201 
202  // These avoid compiler errors; they are only used in sensible context but
203  // compilers typically cannot detect when we access something like
204  // gradients_quad[2] only for dim==3.
205  const unsigned int d1 = dim>1?1:0;
206  const unsigned int d2 = dim>2?2:0;
207  const unsigned int d3 = dim>2?3:0;
208  const unsigned int d4 = dim>2?4:0;
209  const unsigned int d5 = dim>2?5:0;
210 
211  switch (dim)
212  {
213  case 1:
214  for (unsigned int c=0; c<n_components; c++)
215  {
216  if (evaluate_val == true)
217  eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
218  if (evaluate_grad == true)
219  eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
220  if (evaluate_lapl == true)
221  eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
222  }
223  break;
224 
225  case 2:
226  for (unsigned int c=0; c<n_components; c++)
227  {
228  // grad x
229  if (evaluate_grad == true)
230  {
231  eval.template gradients<0,true,false> (values_dofs[c], temp1);
232  eval.template values<1,true,false> (temp1, gradients_quad[c][0]);
233  }
234  if (evaluate_lapl == true)
235  {
236  // grad xy
237  if (evaluate_grad == false)
238  eval.template gradients<0,true,false>(values_dofs[c], temp1);
239  eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]);
240 
241  // grad xx
242  eval.template hessians<0,true,false>(values_dofs[c], temp1);
243  eval.template values<1,true,false> (temp1, hessians_quad[c][0]);
244  }
245 
246  // grad y
247  eval.template values<0,true,false> (values_dofs[c], temp1);
248  if (evaluate_grad == true)
249  eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
250 
251  // grad yy
252  if (evaluate_lapl == true)
253  eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
254 
255  // val: can use values applied in x
256  if (evaluate_val == true)
257  eval.template values<1,true,false> (temp1, values_quad[c]);
258  }
259  break;
260 
261  case 3:
262  for (unsigned int c=0; c<n_components; c++)
263  {
264  if (evaluate_grad == true)
265  {
266  // grad x
267  eval.template gradients<0,true,false> (values_dofs[c], temp1);
268  eval.template values<1,true,false> (temp1, temp2);
269  eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
270  }
271 
272  if (evaluate_lapl == true)
273  {
274  // grad xz
275  if (evaluate_grad == false)
276  {
277  eval.template gradients<0,true,false> (values_dofs[c], temp1);
278  eval.template values<1,true,false> (temp1, temp2);
279  }
280  eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]);
281 
282  // grad xy
283  eval.template gradients<1,true,false> (temp1, temp2);
284  eval.template values<2,true,false> (temp2, hessians_quad[c][d3]);
285 
286  // grad xx
287  eval.template hessians<0,true,false>(values_dofs[c], temp1);
288  eval.template values<1,true,false> (temp1, temp2);
289  eval.template values<2,true,false> (temp2, hessians_quad[c][0]);
290  }
291 
292  // grad y
293  eval.template values<0,true,false> (values_dofs[c], temp1);
294  if (evaluate_grad == true)
295  {
296  eval.template gradients<1,true,false>(temp1, temp2);
297  eval.template values<2,true,false> (temp2, gradients_quad[c][d1]);
298  }
299 
300  if (evaluate_lapl == true)
301  {
302  // grad yz
303  if (evaluate_grad == false)
304  eval.template gradients<1,true,false>(temp1, temp2);
305  eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]);
306 
307  // grad yy
308  eval.template hessians<1,true,false> (temp1, temp2);
309  eval.template values<2,true,false> (temp2, hessians_quad[c][d1]);
310  }
311 
312  // grad z: can use the values applied in x direction stored in temp1
313  eval.template values<1,true,false> (temp1, temp2);
314  if (evaluate_grad == true)
315  eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]);
316 
317  // grad zz: can use the values applied in x and y direction stored
318  // in temp2
319  if (evaluate_lapl == true)
320  eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
321 
322  // val: can use the values applied in x & y direction stored in temp2
323  if (evaluate_val == true)
324  eval.template values<2,true,false> (temp2, values_quad[c]);
325  }
326  break;
327 
328  default:
329  AssertThrow(false, ExcNotImplemented());
330  }
331 
332  // case additional dof for FE_Q_DG0: add values; gradients and second
333  // derivatives evaluate to zero
334  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val)
335  for (unsigned int c=0; c<n_components; ++c)
336  for (unsigned int q=0; q<shape_info.n_q_points; ++q)
337  values_quad[c][q] += values_dofs[c][shape_info.dofs_per_cell-1];
338  }
339 
340 
341 
342  template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
343  int n_q_points_1d, int n_components, typename Number>
344  inline
345  void
346  FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
347  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
348  VectorizedArray<Number> *values_dofs_actual[],
349  VectorizedArray<Number> *values_quad[],
350  VectorizedArray<Number> *gradients_quad[][dim],
351  VectorizedArray<Number> *scratch_data,
352  const bool integrate_val,
353  const bool integrate_grad)
354  {
355  const EvaluatorVariant variant =
356  EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
357  typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
359  Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
360  shape_info.shape_values,
361  variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
362  shape_info.shape_gradients,
363  variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
364  shape_info.shape_hessians,
365  shape_info.fe_degree,
366  shape_info.n_q_points_1d);
367 
368  const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
369  : (Eval::dofs_per_cell > Eval::n_q_points ?
370  Eval::dofs_per_cell : Eval::n_q_points);
371  const unsigned int max_stack_size = 100;
372  VectorizedArray<Number> temp_data[(temp_size > 0 && temp_size < max_stack_size) ? 2*temp_size : 1];
375  if (temp_size == 0)
376  {
377  temp1 = scratch_data;
378  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
379  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
380  }
381  else if (temp_size < max_stack_size)
382  {
383  temp1 = &temp_data[0];
384  temp2 = temp1 + temp_size;
385  }
386  else
387  {
388  temp1 = scratch_data;
389  temp2 = temp1 + temp_size;
390  }
391 
392  // expand dof_values to tensor product for truncated tensor products
393  VectorizedArray<Number> **values_dofs = values_dofs_actual;
394  VectorizedArray<Number> *expanded_dof_values[n_components];
395  if (type == MatrixFreeFunctions::truncated_tensor)
396  {
397  values_dofs = expanded_dof_values;
398  for (unsigned int c=0; c<n_components; ++c)
399  expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
400  shape_info.n_q_points)) +
401  c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
402  }
403 
404  // These avoid compiler errors; they are only used in sensible context but
405  // compilers typically cannot detect when we access something like
406  // gradients_quad[2] only for dim==3.
407  const unsigned int d1 = dim>1?1:0;
408  const unsigned int d2 = dim>2?2:0;
409 
410  switch (dim)
411  {
412  case 1:
413  for (unsigned int c=0; c<n_components; c++)
414  {
415  if (integrate_val == true)
416  eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
417  if (integrate_grad == true)
418  {
419  if (integrate_val == true)
420  eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]);
421  else
422  eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]);
423  }
424  }
425  break;
426 
427  case 2:
428  for (unsigned int c=0; c<n_components; c++)
429  {
430  if (integrate_val == true)
431  {
432  // val
433  eval.template values<0,false,false> (values_quad[c], temp1);
434  //grad x
435  if (integrate_grad == true)
436  eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
437  eval.template values<1,false,false>(temp1, values_dofs[c]);
438  }
439  if (integrate_grad == true)
440  {
441  // grad y
442  eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
443  if (integrate_val == false)
444  {
445  eval.template gradients<1,false,false>(temp1, values_dofs[c]);
446  //grad x
447  eval.template gradients<0,false,false> (gradients_quad[c][0], temp1);
448  eval.template values<1,false,true> (temp1, values_dofs[c]);
449  }
450  else
451  eval.template gradients<1,false,true>(temp1, values_dofs[c]);
452  }
453  }
454  break;
455 
456  case 3:
457  for (unsigned int c=0; c<n_components; c++)
458  {
459  if (integrate_val == true)
460  {
461  // val
462  eval.template values<0,false,false> (values_quad[c], temp1);
463  //grad x: can sum to temporary value in temp1
464  if (integrate_grad == true)
465  eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
466  eval.template values<1,false,false>(temp1, temp2);
467  if (integrate_grad == true)
468  {
469  eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
470  eval.template gradients<1,false,true>(temp1, temp2);
471  }
472  eval.template values<2,false,false> (temp2, values_dofs[c]);
473  }
474  else if (integrate_grad == true)
475  {
476  eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
477  eval.template values<1,false,false> (temp1, temp2);
478  eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
479  eval.template gradients<1,false,true>(temp1, temp2);
480  eval.template values<2,false,false> (temp2, values_dofs[c]);
481  }
482  if (integrate_grad == true)
483  {
484  // grad z: can sum to temporary x and y value in output
485  eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
486  eval.template values<1,false,false> (temp1, temp2);
487  eval.template gradients<2,false,true> (temp2, values_dofs[c]);
488  }
489  }
490  break;
491 
492  default:
493  AssertThrow(false, ExcNotImplemented());
494  }
495 
496  // case FE_Q_DG0: add values, gradients and second derivatives are zero
497  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
498  {
499  if (integrate_val)
500  for (unsigned int c=0; c<n_components; ++c)
501  {
502  values_dofs[c][shape_info.dofs_per_cell-1] = values_quad[c][0];
503  for (unsigned int q=1; q<shape_info.n_q_points; ++q)
504  values_dofs[c][shape_info.dofs_per_cell-1] += values_quad[c][q];
505  }
506  else
507  for (unsigned int c=0; c<n_components; ++c)
508  values_dofs[c][shape_info.dofs_per_cell-1] = VectorizedArray<Number>();
509  }
510 
511  if (type == MatrixFreeFunctions::truncated_tensor)
512  {
513  unsigned int count_p = 0, count_q = 0;
514  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
515  for (int i=0; i<(dim>2?degree+1:1); ++i)
516  {
517  for (int j=0; j<(dim>1?degree+1-i:1); ++j)
518  {
519  for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
520  {
521  for (unsigned int c=0; c<n_components; ++c)
522  values_dofs_actual[c][count_p] = expanded_dof_values[c][count_q];
523  }
524  count_q += j+i;
525  }
526  count_q += i*(degree+1);
527  }
528  AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
529  }
530  }
531 
532  // This a specialization for Gauss-Lobatto elements where the 'values'
533  // operation is identity, which allows us to write shorter code.
534  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
535  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
536  fe_degree, n_q_points_1d, n_components, Number>
537  {
538  static
539  void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
540  VectorizedArray<Number> *values_dofs[],
541  VectorizedArray<Number> *values_quad[],
542  VectorizedArray<Number> *gradients_quad[][dim],
543  VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
544  VectorizedArray<Number> *scratch_data,
545  const bool evaluate_val,
546  const bool evaluate_grad,
547  const bool evaluate_lapl);
548 
549  static
550  void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
551  VectorizedArray<Number> *values_dofs[],
552  VectorizedArray<Number> *values_quad[],
553  VectorizedArray<Number> *gradients_quad[][dim],
554  VectorizedArray<Number> *scratch_data,
555  const bool integrate_val,
556  const bool integrate_grad);
557  };
558 
559  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
560  inline
561  void
562  FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
563  fe_degree, n_q_points_1d, n_components, Number>
564  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
565  VectorizedArray<Number> *values_dofs[],
566  VectorizedArray<Number> *values_quad[],
567  VectorizedArray<Number> *gradients_quad[][dim],
568  VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
569  VectorizedArray<Number> *scratch_data,
570  const bool evaluate_val,
571  const bool evaluate_grad,
572  const bool evaluate_lapl)
573  {
574  typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
576  Eval eval (shape_info.shape_val_evenodd,
577  shape_info.shape_gra_evenodd,
578  shape_info.shape_hes_evenodd,
579  shape_info.fe_degree,
580  shape_info.n_q_points_1d);
581 
582  // These avoid compiler errors; they are only used in sensible context but
583  // compilers typically cannot detect when we access something like
584  // gradients_quad[2] only for dim==3.
585  const unsigned int d1 = dim>1?1:0;
586  const unsigned int d2 = dim>2?2:0;
587  const unsigned int d3 = dim>2?3:0;
588  const unsigned int d4 = dim>2?4:0;
589  const unsigned int d5 = dim>2?5:0;
590 
591  switch (dim)
592  {
593  case 1:
594  if (evaluate_val == true)
595  std::memcpy (values_quad[0], values_dofs[0],
596  eval.dofs_per_cell * n_components *
597  sizeof (values_dofs[0][0]));
598  for (unsigned int c=0; c<n_components; c++)
599  {
600  if (evaluate_grad == true)
601  eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
602  if (evaluate_lapl == true)
603  eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
604  }
605  break;
606 
607  case 2:
608  if (evaluate_val == true)
609  {
610  std::memcpy (values_quad[0], values_dofs[0],
611  Eval::dofs_per_cell * n_components *
612  sizeof (values_dofs[0][0]));
613  }
614  if (evaluate_grad == true)
615  for (unsigned int comp=0; comp<n_components; comp++)
616  {
617  // grad x
618  eval.template gradients<0,true,false> (values_dofs[comp],
619  gradients_quad[comp][0]);
620  // grad y
621  eval.template gradients<1,true,false> (values_dofs[comp],
622  gradients_quad[comp][d1]);
623  }
624  if (evaluate_lapl == true)
625  for (unsigned int comp=0; comp<n_components; comp++)
626  {
627  // hess x
628  eval.template hessians<0,true,false> (values_dofs[comp],
629  hessians_quad[comp][0]);
630  // hess y
631  eval.template hessians<1,true,false> (values_dofs[comp],
632  hessians_quad[comp][d1]);
633 
634  // grad x grad y
635  eval.template gradients<0,true,false> (values_dofs[comp], scratch_data);
636  eval.template gradients<1,true,false> (scratch_data, hessians_quad[comp][d1+d1]);
637  }
638  break;
639 
640  case 3:
641  if (evaluate_val == true)
642  {
643  std::memcpy (values_quad[0], values_dofs[0],
644  Eval::dofs_per_cell * n_components *
645  sizeof (values_dofs[0][0]));
646  }
647  if (evaluate_grad == true)
648  for (unsigned int comp=0; comp<n_components; comp++)
649  {
650  // grad x
651  eval.template gradients<0,true,false> (values_dofs[comp],
652  gradients_quad[comp][0]);
653  // grad y
654  eval.template gradients<1,true,false> (values_dofs[comp],
655  gradients_quad[comp][d1]);
656  // grad y
657  eval.template gradients<2,true,false> (values_dofs[comp],
658  gradients_quad[comp][d2]);
659  }
660  if (evaluate_lapl == true)
661  for (unsigned int comp=0; comp<n_components; comp++)
662  {
663  // grad x
664  eval.template hessians<0,true,false> (values_dofs[comp],
665  hessians_quad[comp][0]);
666  // grad y
667  eval.template hessians<1,true,false> (values_dofs[comp],
668  hessians_quad[comp][d1]);
669  // grad y
670  eval.template hessians<2,true,false> (values_dofs[comp],
671  hessians_quad[comp][d2]);
672 
673  VectorizedArray<Number> *temp1 = scratch_data;
674  // grad xy
675  eval.template gradients<0,true,false> (values_dofs[comp], temp1);
676  eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]);
677  // grad xz
678  eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]);
679  // grad yz
680  eval.template gradients<1,true,false> (values_dofs[comp], temp1);
681  eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]);
682  }
683  break;
684  default:
685  AssertThrow(false, ExcNotImplemented());
686  }
687  }
688 
689  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
690  inline
691  void
692  FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
693  fe_degree, n_q_points_1d, n_components, Number>
694  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
695  VectorizedArray<Number> *values_dofs[],
696  VectorizedArray<Number> *values_quad[],
697  VectorizedArray<Number> *gradients_quad[][dim],
699  const bool integrate_val,
700  const bool integrate_grad)
701  {
702  typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
704  Eval eval (shape_info.shape_val_evenodd,
705  shape_info.shape_gra_evenodd,
706  shape_info.shape_hes_evenodd,
707  shape_info.fe_degree,
708  shape_info.n_q_points_1d);
709 
710  // These avoid compiler errors; they are only used in sensible context but
711  // compilers typically cannot detect when we access something like
712  // gradients_quad[2] only for dim==3.
713  const unsigned int d1 = dim>1?1:0;
714  const unsigned int d2 = dim>2?2:0;
715 
716  if (integrate_val == true)
717  std::memcpy (values_dofs[0], values_quad[0],
718  Eval::dofs_per_cell * n_components *
719  sizeof (values_dofs[0][0]));
720  switch (dim)
721  {
722  case 1:
723  for (unsigned int c=0; c<n_components; c++)
724  {
725  if (integrate_grad == true)
726  {
727  if (integrate_val == true)
728  eval.template gradients<0,false,true> (gradients_quad[c][0],
729  values_dofs[c]);
730  else
731  eval.template gradients<0,false,false> (gradients_quad[c][0],
732  values_dofs[c]);
733  }
734  }
735 
736  break;
737  case 2:
738  if (integrate_grad == true)
739  for (unsigned int comp=0; comp<n_components; comp++)
740  {
741  // grad x: If integrate_val == true we have to add to the
742  // previous output
743  if (integrate_val == true)
744  eval.template gradients<0, false, true> (gradients_quad[comp][0],
745  values_dofs[comp]);
746  else
747  eval.template gradients<0, false, false> (gradients_quad[comp][0],
748  values_dofs[comp]);
749 
750  // grad y
751  eval.template gradients<1, false, true> (gradients_quad[comp][d1],
752  values_dofs[comp]);
753  }
754  break;
755 
756  case 3:
757  if (integrate_grad == true)
758  for (unsigned int comp=0; comp<n_components; comp++)
759  {
760  // grad x: If integrate_val == true we have to add to the
761  // previous output
762  if (integrate_val == true)
763  eval.template gradients<0, false, true> (gradients_quad[comp][0],
764  values_dofs[comp]);
765  else
766  eval.template gradients<0, false, false> (gradients_quad[comp][0],
767  values_dofs[comp]);
768 
769  // grad y
770  eval.template gradients<1, false, true> (gradients_quad[comp][d1],
771  values_dofs[comp]);
772 
773  // grad z
774  eval.template gradients<2, false, true> (gradients_quad[comp][d2],
775  values_dofs[comp]);
776  }
777  break;
778 
779  default:
780  AssertThrow(false, ExcNotImplemented());
781  }
782  }
783 
784 } // end of namespace internal
785 
786 
787 DEAL_II_NAMESPACE_CLOSE
788 
789 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:170
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1146
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
static ::ExceptionBase & ExcNotImplemented()