Reference documentation for deal.II version 9.0.0
tensor_product_kernels.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018-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_tensor_product_kernels_h
18 #define dealii_matrix_free_tensor_product_kernels_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/aligned_vector.h>
22 #include <deal.II/base/utilities.h>
23 
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 
29 namespace internal
30 {
35  enum EvaluatorVariant
36  {
41  evaluate_general,
47  evaluate_symmetric,
53  evaluate_evenodd,
64  evaluate_symmetric_hierarchical
65  };
66 
67 
68 
89  template <EvaluatorVariant variant, int dim, int n_rows, int n_columns,
90  typename Number, typename Number2=Number>
91  struct EvaluatorTensorProduct
92  {};
93 
94 
95 
113  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
114  struct EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>
115  {
116  static constexpr unsigned int n_rows_of_product = Utilities::pow(n_rows, dim);
117  static constexpr unsigned int n_columns_of_product = Utilities::pow(n_columns, dim);
118 
123  EvaluatorTensorProduct ()
124  :
125  shape_values (nullptr),
126  shape_gradients (nullptr),
127  shape_hessians (nullptr)
128  {}
129 
133  EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
134  const AlignedVector<Number2> &shape_gradients,
135  const AlignedVector<Number2> &shape_hessians,
136  const unsigned int dummy1 = 0,
137  const unsigned int dummy2 = 0)
138  :
139  shape_values (shape_values.begin()),
140  shape_gradients (shape_gradients.begin()),
141  shape_hessians (shape_hessians.begin())
142  {
143  // We can enter this function either for the apply() path that has
144  // n_rows * n_columns entries or for the apply_face() path that only has
145  // n_rows * 3 entries in the array. Since we cannot decide about the use
146  // we must allow for both here.
147  Assert(shape_values.size() == 0 ||
148  shape_values.size() == n_rows*n_columns ||
149  shape_values.size() == 3*n_rows,
150  ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));
151  Assert(shape_gradients.size() == 0 ||
152  shape_gradients.size() == n_rows*n_columns,
153  ExcDimensionMismatch(shape_gradients.size(), n_rows*n_columns));
154  Assert(shape_hessians.size() == 0 ||
155  shape_hessians.size() == n_rows*n_columns,
156  ExcDimensionMismatch(shape_hessians.size(), n_rows*n_columns));
157  (void)dummy1;
158  (void)dummy2;
159  }
160 
161  template <int direction, bool contract_over_rows, bool add>
162  void
163  values (const Number in [],
164  Number out[]) const
165  {
166  apply<direction,contract_over_rows,add>(shape_values, in, out);
167  }
168 
169  template <int direction, bool contract_over_rows, bool add>
170  void
171  gradients (const Number in [],
172  Number out[]) const
173  {
174  apply<direction,contract_over_rows,add>(shape_gradients, in, out);
175  }
176 
177  template <int direction, bool contract_over_rows, bool add>
178  void
179  hessians (const Number in [],
180  Number out[]) const
181  {
182  apply<direction,contract_over_rows,add>(shape_hessians, in, out);
183  }
184 
185  template <int direction, bool contract_over_rows, bool add>
186  void
187  values_one_line (const Number in [],
188  Number out[]) const
189  {
190  Assert (shape_values != nullptr, ExcNotInitialized());
191  apply<direction,contract_over_rows,add,true>(shape_values, in, out);
192  }
193 
194  template <int direction, bool contract_over_rows, bool add>
195  void
196  gradients_one_line (const Number in [],
197  Number out[]) const
198  {
199  Assert (shape_gradients != nullptr, ExcNotInitialized());
200  apply<direction,contract_over_rows,add,true>(shape_gradients, in, out);
201  }
202 
203  template <int direction, bool contract_over_rows, bool add>
204  void
205  hessians_one_line (const Number in [],
206  Number out[]) const
207  {
208  Assert (shape_hessians != nullptr, ExcNotInitialized());
209  apply<direction,contract_over_rows,add,true>(shape_hessians, in, out);
210  }
211 
236  template <int direction, bool contract_over_rows, bool add, bool one_line=false>
237  static void apply (const Number2 *DEAL_II_RESTRICT shape_data,
238  const Number *in,
239  Number *out);
240 
270  template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
271  void apply_face (const Number *DEAL_II_RESTRICT in,
272  Number *DEAL_II_RESTRICT out) const;
273 
274  const Number2 *shape_values;
275  const Number2 *shape_gradients;
276  const Number2 *shape_hessians;
277  };
278 
279 
280 
281  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
282  template <int direction, bool contract_over_rows, bool add, bool one_line>
283  inline
284  void
285  EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>
286  ::apply (const Number2 *DEAL_II_RESTRICT shape_data,
287  const Number *in,
288  Number *out)
289  {
290  static_assert (one_line == false || direction==dim-1,
291  "Single-line evaluation only works for direction=dim-1.");
292  Assert(shape_data != nullptr,
293  ExcMessage("The given array shape_data must not be the null pointer!"));
294  Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
295  ExcMessage("In-place operation only supported for "
296  "n_rows==n_columns or single-line interpolation"));
297  AssertIndexRange (direction, dim);
298  constexpr int mm = contract_over_rows ? n_rows : n_columns,
299  nn = contract_over_rows ? n_columns : n_rows;
300 
301  constexpr int stride = Utilities::pow(n_columns, direction);
302  constexpr int n_blocks1 = one_line ? 1 : stride;
303  constexpr int n_blocks2 = Utilities::pow(n_rows, (direction>=dim) ? 0 : (dim-direction-1));
304 
305  for (int i2=0; i2<n_blocks2; ++i2)
306  {
307  for (int i1=0; i1<n_blocks1; ++i1)
308  {
309  Number x[mm];
310  for (int i=0; i<mm; ++i)
311  x[i] = in[stride*i];
312  for (int col=0; col<nn; ++col)
313  {
314  Number2 val0;
315  if (contract_over_rows == true)
316  val0 = shape_data[col];
317  else
318  val0 = shape_data[col*n_columns];
319  Number res0 = val0 * x[0];
320  for (int i=1; i<mm; ++i)
321  {
322  if (contract_over_rows == true)
323  val0 = shape_data[i*n_columns+col];
324  else
325  val0 = shape_data[col*n_columns+i];
326  res0 += val0 * x[i];
327  }
328  if (add == false)
329  out[stride*col] = res0;
330  else
331  out[stride*col] += res0;
332  }
333 
334  if (one_line == false)
335  {
336  ++in;
337  ++out;
338  }
339  }
340  if (one_line == false)
341  {
342  in += stride*(mm-1);
343  out += stride*(nn-1);
344  }
345  }
346  }
347 
348 
349 
350  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
351  template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
352  inline
353  void
354  EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>
355  ::apply_face (const Number *DEAL_II_RESTRICT in,
356  Number *DEAL_II_RESTRICT out) const
357  {
358  static_assert(dim > 0 && dim<4, "Only dim=1,2,3 supported");
359  static_assert(max_derivative >= 0 && max_derivative<3,
360  "Only derivative orders 0-2 implemented");
361  Assert(shape_values != nullptr,
362  ExcMessage("The given array shape_values must not be the null pointer."));
363 
364  constexpr int n_blocks1 = dim > 1 ? n_rows : 1;
365  constexpr int n_blocks2 = dim > 2 ? n_rows : 1;
366 
367  AssertIndexRange (face_direction, dim);
368  constexpr int stride = Utilities::pow(n_rows, face_direction);
369  constexpr int out_stride = Utilities::pow(n_rows, dim-1);
370  const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
371 
372  for (int i2=0; i2<n_blocks2; ++i2)
373  {
374  for (int i1=0; i1<n_blocks1; ++i1)
375  {
376  if (contract_onto_face == true)
377  {
378  Number res0 = shape_values[0] * in[0];
379  Number res1, res2;
380  if (max_derivative > 0)
381  res1 = shape_values[n_rows] * in[0];
382  if (max_derivative > 1)
383  res2 = shape_values[2*n_rows] * in[0];
384  for (int ind=1; ind<n_rows; ++ind)
385  {
386  res0 += shape_values[ind] * in[stride*ind];
387  if (max_derivative > 0)
388  res1 += shape_values[ind+n_rows] * in[stride*ind];
389  if (max_derivative > 1)
390  res2 += shape_values[ind+2*n_rows] * in[stride*ind];
391  }
392  if (add == false)
393  {
394  out[0] = res0;
395  if (max_derivative > 0)
396  out[out_stride] = res1;
397  if (max_derivative > 1)
398  out[2*out_stride] = res2;
399  }
400  else
401  {
402  out[0] += res0;
403  if (max_derivative > 0)
404  out[out_stride] += res1;
405  if (max_derivative > 1)
406  out[2*out_stride] += res2;
407  }
408  }
409  else
410  {
411  for (int col=0; col<n_rows; ++col)
412  {
413  if (add == false)
414  out[col*stride] = shape_values[col] * in[0];
415  else
416  out[col*stride] += shape_values[col] * in[0];
417  if (max_derivative > 0)
418  out[col*stride] += shape_values[col+n_rows] * in[out_stride];
419  if (max_derivative > 1)
420  out[col*stride] += shape_values[col+2*n_rows] * in[2*out_stride];
421  }
422  }
423 
424  // increment: in regular case, just go to the next point in
425  // x-direction. If we are at the end of one chunk in x-dir, need
426  // to jump over to the next layer in z-direction
427  switch (face_direction)
428  {
429  case 0:
430  in += contract_onto_face ? n_rows : 1;
431  out += contract_onto_face ? 1 : n_rows;
432  break;
433  case 1:
434  ++in;
435  ++out;
436  // faces 2 and 3 in 3D use local coordinate system zx, which
437  // is the other way around compared to the tensor
438  // product. Need to take that into account.
439  if (dim == 3)
440  {
441  if (contract_onto_face)
442  out += n_rows-1;
443  else
444  in += n_rows-1;
445  }
446  break;
447  case 2:
448  ++in;
449  ++out;
450  break;
451  default:
452  Assert (false, ExcNotImplemented());
453  }
454  }
455  if (face_direction == 1 && dim == 3)
456  {
457  // adjust for local coordinate system zx
458  if (contract_onto_face)
459  {
460  in += n_rows*(n_rows-1);
461  out -= n_rows*n_rows-1;
462  }
463  else
464  {
465  out += n_rows*(n_rows-1);
466  in -= n_rows*n_rows-1;
467  }
468  }
469  }
470  }
471 
472 
473 
487  template <int dim, typename Number, typename Number2>
488  struct EvaluatorTensorProduct<evaluate_general,dim,0,0,Number,Number2>
489  {
490  static constexpr unsigned int n_rows_of_product = numbers::invalid_unsigned_int;
491  static constexpr unsigned int n_columns_of_product = numbers::invalid_unsigned_int;
492 
497  EvaluatorTensorProduct ()
498  :
499  shape_values (nullptr),
500  shape_gradients (nullptr),
501  shape_hessians (nullptr),
502  n_rows (numbers::invalid_unsigned_int),
503  n_columns (numbers::invalid_unsigned_int)
504  {}
505 
509  EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
510  const AlignedVector<Number2> &shape_gradients,
511  const AlignedVector<Number2> &shape_hessians,
512  const unsigned int n_rows,
513  const unsigned int n_columns)
514  :
515  shape_values (shape_values.begin()),
516  shape_gradients (shape_gradients.begin()),
517  shape_hessians (shape_hessians.begin()),
518  n_rows (n_rows),
519  n_columns (n_columns)
520  {
521  // We can enter this function either for the apply() path that has
522  // n_rows * n_columns entries or for the apply_face() path that only has
523  // n_rows * 3 entries in the array. Since we cannot decide about the use
524  // we must allow for both here.
525  Assert(shape_values.size() == 0 ||
526  shape_values.size() == n_rows*n_columns ||
527  shape_values.size() == n_rows*3,
528  ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));
529  Assert(shape_gradients.size() == 0 ||
530  shape_gradients.size() == n_rows*n_columns,
531  ExcDimensionMismatch(shape_gradients.size(), n_rows*n_columns));
532  Assert(shape_hessians.size() == 0 ||
533  shape_hessians.size() == n_rows*n_columns,
534  ExcDimensionMismatch(shape_hessians.size(), n_rows*n_columns));
535  }
536 
537  template <int direction, bool contract_over_rows, bool add>
538  void
539  values (const Number *in,
540  Number *out) const
541  {
542  apply<direction,contract_over_rows,add>(shape_values, in, out);
543  }
544 
545  template <int direction, bool contract_over_rows, bool add>
546  void
547  gradients (const Number *in,
548  Number *out) const
549  {
550  apply<direction,contract_over_rows,add>(shape_gradients, in, out);
551  }
552 
553  template <int direction, bool contract_over_rows, bool add>
554  void
555  hessians (const Number *in,
556  Number *out) const
557  {
558  apply<direction,contract_over_rows,add>(shape_hessians, in, out);
559  }
560 
561  template <int direction, bool contract_over_rows, bool add>
562  void
563  values_one_line (const Number in [],
564  Number out[]) const
565  {
566  Assert (shape_values != nullptr, ExcNotInitialized());
567  apply<direction,contract_over_rows,add,true>(shape_values, in, out);
568  }
569 
570  template <int direction, bool contract_over_rows, bool add>
571  void
572  gradients_one_line (const Number in [],
573  Number out[]) const
574  {
575  Assert (shape_gradients != nullptr, ExcNotInitialized());
576  apply<direction,contract_over_rows,add,true>(shape_gradients, in, out);
577  }
578 
579  template <int direction, bool contract_over_rows, bool add>
580  void
581  hessians_one_line (const Number in [],
582  Number out[]) const
583  {
584  Assert (shape_hessians != nullptr, ExcNotInitialized());
585  apply<direction,contract_over_rows,add,true>(shape_hessians, in, out);
586  }
587 
588  template <int direction, bool contract_over_rows, bool add, bool one_line=false>
589  void apply (const Number2 *DEAL_II_RESTRICT shape_data,
590  const Number *in,
591  Number *out) const;
592 
593  template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
594  void apply_face (const Number *DEAL_II_RESTRICT in,
595  Number *DEAL_II_RESTRICT out) const;
596 
597  const Number2 *shape_values;
598  const Number2 *shape_gradients;
599  const Number2 *shape_hessians;
600  const unsigned int n_rows;
601  const unsigned int n_columns;
602  };
603 
604 
605 
606  template <int dim, typename Number, typename Number2>
607  template <int direction, bool contract_over_rows, bool add, bool one_line>
608  inline
609  void
610  EvaluatorTensorProduct<evaluate_general,dim,0,0,Number,Number2>
611  ::apply (const Number2 *DEAL_II_RESTRICT shape_data,
612  const Number *in,
613  Number *out) const
614  {
615  static_assert (one_line == false || direction==dim-1,
616  "Single-line evaluation only works for direction=dim-1.");
617  Assert(shape_data != nullptr,
618  ExcMessage("The given array shape_data must not be the null pointer!"));
619  Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
620  ExcMessage("In-place operation only supported for "
621  "n_rows==n_columns or single-line interpolation"));
622  AssertIndexRange (direction, dim);
623  const int mm = contract_over_rows ? n_rows : n_columns,
624  nn = contract_over_rows ? n_columns : n_rows;
625 
626  const int stride = direction==0 ? 1 : Utilities::fixed_power<direction>(n_columns);
627  const int n_blocks1 = one_line ? 1 : stride;
628  const int n_blocks2 = direction >= dim-1 ? 1 : Utilities::fixed_power<dim-direction-1>(n_rows);
629  Assert(n_rows <= 128, ExcNotImplemented());
630 
631  for (int i2=0; i2<n_blocks2; ++i2)
632  {
633  for (int i1=0; i1<n_blocks1; ++i1)
634  {
635  Number x[129];
636  for (int i=0; i<mm; ++i)
637  x[i] = in[stride*i];
638  for (int col=0; col<nn; ++col)
639  {
640  Number2 val0;
641  if (contract_over_rows == true)
642  val0 = shape_data[col];
643  else
644  val0 = shape_data[col*n_columns];
645  Number res0 = val0 * x[0];
646  for (int i=1; i<mm; ++i)
647  {
648  if (contract_over_rows == true)
649  val0 = shape_data[i*n_columns+col];
650  else
651  val0 = shape_data[col*n_columns+i];
652  res0 += val0 * x[i];
653  }
654  if (add == false)
655  out[stride*col] = res0;
656  else
657  out[stride*col] += res0;
658  }
659 
660  if (one_line == false)
661  {
662  ++in;
663  ++out;
664  }
665  }
666  if (one_line == false)
667  {
668  in += stride*(mm-1);
669  out += stride*(nn-1);
670  }
671  }
672  }
673 
674 
675 
676  template <int dim, typename Number, typename Number2>
677  template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
678  inline
679  void
680  EvaluatorTensorProduct<evaluate_general,dim,0,0,Number,Number2>
681  ::apply_face (const Number *DEAL_II_RESTRICT in,
682  Number *DEAL_II_RESTRICT out) const
683  {
684  Assert(shape_values != nullptr,
685  ExcMessage("The given array shape_data must not be the null pointer!"));
686  static_assert(dim > 0 && dim<4, "Only dim=1,2,3 supported");
687  const int n_blocks1 = dim > 1 ? n_rows : 1;
688  const int n_blocks2 = dim > 2 ? n_rows : 1;
689 
690  AssertIndexRange (face_direction, dim);
691  const int stride = face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
692  const int out_stride = dim > 1 ? Utilities::fixed_power<dim-1>(n_rows) : 1;
693 
694  for (int i2=0; i2<n_blocks2; ++i2)
695  {
696  for (int i1=0; i1<n_blocks1; ++i1)
697  {
698  if (contract_onto_face == true)
699  {
700  Number res0 = shape_values[0] * in[0];
701  Number res1, res2;
702  if (max_derivative > 0)
703  res1 = shape_values[n_rows] * in[0];
704  if (max_derivative > 1)
705  res2 = shape_values[2*n_rows] * in[0];
706  for (unsigned int ind=1; ind<n_rows; ++ind)
707  {
708  res0 += shape_values[ind] * in[stride*ind];
709  if (max_derivative > 0)
710  res1 += shape_values[ind+n_rows] * in[stride*ind];
711  if (max_derivative > 1)
712  res2 += shape_values[ind+2*n_rows] * in[stride*ind];
713  }
714  if (add == false)
715  {
716  out[0] = res0;
717  if (max_derivative > 0)
718  out[out_stride] = res1;
719  if (max_derivative > 1)
720  out[2*out_stride] = res2;
721  }
722  else
723  {
724  out[0] += res0;
725  if (max_derivative > 0)
726  out[out_stride] += res1;
727  if (max_derivative > 1)
728  out[2*out_stride] += res2;
729  }
730  }
731  else
732  {
733  for (unsigned int col=0; col<n_rows; ++col)
734  {
735  if (add == false)
736  out[col*stride] = shape_values[col] * in[0];
737  else
738  out[col*stride] += shape_values[col] * in[0];
739  if (max_derivative > 0)
740  out[col*stride] += shape_values[col+n_rows] * in[out_stride];
741  if (max_derivative > 1)
742  out[col*stride] += shape_values[col+2*n_rows] * in[2*out_stride];
743  }
744  }
745 
746  // increment: in regular case, just go to the next point in
747  // x-direction. If we are at the end of one chunk in x-dir, need
748  // to jump over to the next layer in z-direction
749  switch (face_direction)
750  {
751  case 0:
752  in += contract_onto_face ? n_rows : 1;
753  out += contract_onto_face ? 1 : n_rows;
754  break;
755  case 1:
756  ++in;
757  ++out;
758  // faces 2 and 3 in 3D use local coordinate system zx, which
759  // is the other way around compared to the tensor
760  // product. Need to take that into account.
761  if (dim == 3)
762  {
763  if (contract_onto_face)
764  out += n_rows-1;
765  else
766  in += n_rows-1;
767  }
768  break;
769  case 2:
770  ++in;
771  ++out;
772  break;
773  default:
774  Assert (false, ExcNotImplemented());
775  }
776  }
777  if (face_direction == 1 && dim == 3)
778  {
779  // adjust for local coordinate system zx
780  if (contract_onto_face)
781  {
782  in += n_rows*(n_rows-1);
783  out -= n_rows*n_rows-1;
784  }
785  else
786  {
787  out += n_rows*(n_rows-1);
788  in -= n_rows*n_rows-1;
789  }
790  }
791  }
792  }
793 
794 
795 
816  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
817  struct EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
818  {
819  static constexpr unsigned int n_rows_of_product = Utilities::pow(n_rows, dim);
820  static constexpr unsigned int n_columns_of_product = Utilities::pow(n_columns, dim);
821 
825  EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
826  const AlignedVector<Number2> &shape_gradients,
827  const AlignedVector<Number2> &shape_hessians,
828  const unsigned int dummy1 = 0,
829  const unsigned int dummy2 = 0)
830  :
831  shape_values (shape_values.begin()),
832  shape_gradients (shape_gradients.begin()),
833  shape_hessians (shape_hessians.begin())
834  {
835  Assert(shape_values.size() == 0 ||
836  shape_values.size() == n_rows*n_columns,
837  ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));
838  Assert(shape_gradients.size() == 0 ||
839  shape_gradients.size() == n_rows*n_columns,
840  ExcDimensionMismatch(shape_gradients.size(), n_rows*n_columns));
841  Assert(shape_hessians.size() == 0 ||
842  shape_hessians.size() == n_rows*n_columns,
843  ExcDimensionMismatch(shape_hessians.size(), n_rows*n_columns));
844  (void)dummy1;
845  (void)dummy2;
846  }
847 
848  template <int direction, bool contract_over_rows, bool add>
849  void
850  values (const Number in [],
851  Number out[]) const;
852 
853  template <int direction, bool contract_over_rows, bool add>
854  void
855  gradients (const Number in [],
856  Number out[]) const;
857 
858  template <int direction, bool contract_over_rows, bool add>
859  void
860  hessians (const Number in [],
861  Number out[]) const;
862 
863  const Number2 *shape_values;
864  const Number2 *shape_gradients;
865  const Number2 *shape_hessians;
866  };
867 
868 
869 
870  // In this case, the 1D shape values read (sorted lexicographically, rows
871  // run over 1D dofs, columns over quadrature points):
872  // Q2 --> [ 0.687 0 -0.087 ]
873  // [ 0.4 1 0.4 ]
874  // [-0.087 0 0.687 ]
875  // Q3 --> [ 0.66 0.003 0.002 0.049 ]
876  // [ 0.521 1.005 -0.01 -0.230 ]
877  // [-0.230 -0.01 1.005 0.521 ]
878  // [ 0.049 0.002 0.003 0.66 ]
879  // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
880  // [ 0.608 1.059 0 0.039 0.176 ]
881  // [-0.409 -0.113 1 -0.113 -0.409 ]
882  // [ 0.176 0.039 0 1.059 0.608 ]
883  // [-0.032 -0.007 0 0.022 0.658 ]
884  //
885  // In these matrices, we want to use avoid computations involving zeros and
886  // ones and in addition use the symmetry in entries to reduce the number of
887  // read operations.
888  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
889  template <int direction, bool contract_over_rows, bool add>
890  inline
891  void
892  EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
893  ::values (const Number in [],
894  Number out []) const
895  {
896  Assert (shape_values != nullptr, ExcNotInitialized());
897  AssertIndexRange (direction, dim);
898  constexpr int mm = contract_over_rows ? n_rows : n_columns,
899  nn = contract_over_rows ? n_columns : n_rows;
900  constexpr int n_cols = nn / 2;
901  constexpr int mid = mm / 2;
902 
903  constexpr int stride = Utilities::pow(n_columns, direction);
904  constexpr int n_blocks1 = stride;
905  constexpr int n_blocks2 = Utilities::pow(n_rows, (direction>=dim) ? 0 : (dim-direction-1));
906 
907  for (int i2=0; i2<n_blocks2; ++i2)
908  {
909  for (int i1=0; i1<n_blocks1; ++i1)
910  {
911  for (int col=0; col<n_cols; ++col)
912  {
913  Number2 val0, val1;
914  Number in0, in1, res0, res1;
915  if (contract_over_rows == true)
916  {
917  val0 = shape_values[col];
918  val1 = shape_values[nn-1-col];
919  }
920  else
921  {
922  val0 = shape_values[col*n_columns];
923  val1 = shape_values[(col+1)*n_columns-1];
924  }
925  if (mid > 0)
926  {
927  in0 = in[0];
928  in1 = in[stride*(mm-1)];
929  res0 = val0 * in0;
930  res1 = val1 * in0;
931  res0 += val1 * in1;
932  res1 += val0 * in1;
933  for (int ind=1; ind<mid; ++ind)
934  {
935  if (contract_over_rows == true)
936  {
937  val0 = shape_values[ind*n_columns+col];
938  val1 = shape_values[ind*n_columns+nn-1-col];
939  }
940  else
941  {
942  val0 = shape_values[col*n_columns+ind];
943  val1 = shape_values[(col+1)*n_columns-1-ind];
944  }
945  in0 = in[stride*ind];
946  in1 = in[stride*(mm-1-ind)];
947  res0 += val0 * in0;
948  res1 += val1 * in0;
949  res0 += val1 * in1;
950  res1 += val0 * in1;
951  }
952  }
953  else
954  res0 = res1 = Number();
955  if (contract_over_rows == true)
956  {
957  if (mm % 2 == 1)
958  {
959  val0 = shape_values[mid*n_columns+col];
960  in1 = val0 * in[stride*mid];
961  res0 += in1;
962  res1 += in1;
963  }
964  }
965  else
966  {
967  if (mm % 2 == 1 && nn % 2 == 0)
968  {
969  val0 = shape_values[col*n_columns+mid];
970  in1 = val0 * in[stride*mid];
971  res0 += in1;
972  res1 += in1;
973  }
974  }
975  if (add == false)
976  {
977  out[stride*col] = res0;
978  out[stride*(nn-1-col)] = res1;
979  }
980  else
981  {
982  out[stride*col] += res0;
983  out[stride*(nn-1-col)] += res1;
984  }
985  }
986  if ( contract_over_rows == true && nn%2==1 && mm%2==1 )
987  {
988  if (add==false)
989  out[stride*n_cols] = in[stride*mid];
990  else
991  out[stride*n_cols] += in[stride*mid];
992  }
993  else if (contract_over_rows == true && nn%2==1)
994  {
995  Number res0;
996  Number2 val0 = shape_values[n_cols];
997  if (mid > 0)
998  {
999  res0 = val0 * (in[0] + in[stride*(mm-1)]);
1000  for (int ind=1; ind<mid; ++ind)
1001  {
1002  val0 = shape_values[ind*n_columns+n_cols];
1003  res0 += val0 * (in[stride*ind] + in[stride*(mm-1-ind)]);
1004  }
1005  }
1006  else
1007  res0 = Number();
1008  if (add == false)
1009  out[stride*n_cols] = res0;
1010  else
1011  out[stride*n_cols] += res0;
1012  }
1013  else if (contract_over_rows == false && nn%2 == 1)
1014  {
1015  Number res0;
1016  if (mid > 0)
1017  {
1018  Number2 val0 = shape_values[n_cols*n_columns];
1019  res0 = val0 * (in[0] + in[stride*(mm-1)]);
1020  for (int ind=1; ind<mid; ++ind)
1021  {
1022  val0 = shape_values[n_cols*n_columns+ind];
1023  Number in1 = val0 * (in[stride*ind] + in[stride*(mm-1-ind)]);
1024  res0 += in1;
1025  }
1026  if (mm % 2)
1027  res0 += in[stride*mid];
1028  }
1029  else
1030  res0 = in[0];
1031  if (add == false)
1032  out[stride*n_cols] = res0;
1033  else
1034  out[stride*n_cols] += res0;
1035  }
1036 
1037  ++in;
1038  ++out;
1039  }
1040  in += stride*(mm-1);
1041  out += stride*(nn-1);
1042  }
1043  }
1044 
1045 
1046 
1047  // For the specialized loop used for the gradient computation in
1048  // here, the 1D shape values read (sorted lexicographically, rows
1049  // run over 1D dofs, columns over quadrature points):
1050  // Q2 --> [-2.549 -1 0.549 ]
1051  // [ 3.098 0 -3.098 ]
1052  // [-0.549 1 2.549 ]
1053  // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1054  // [ 6.07 -1.44 -2.97 2.196 ]
1055  // [-2.196 2.97 1.44 -6.07 ]
1056  // [ 0.44 -0.5 1.03 4.315 ]
1057  // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1058  // [10.111 -2.76 -2.667 2.066 -2.306 ]
1059  // [-5.688 5.773 0 -5.773 5.688 ]
1060  // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1061  // [-0.413 0.353 -0.333 -0.353 0.413 ]
1062  //
1063  // In these matrices, we want to use avoid computations involving
1064  // zeros and ones and in addition use the symmetry in entries to
1065  // reduce the number of read operations.
1066  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
1067  template <int direction, bool contract_over_rows, bool add>
1068  inline
1069  void
1070  EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
1071  ::gradients (const Number in [],
1072  Number out []) const
1073  {
1074  Assert (shape_gradients != nullptr, ExcNotInitialized());
1075  AssertIndexRange (direction, dim);
1076  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1077  nn = contract_over_rows ? n_columns : n_rows;
1078  constexpr int n_cols = nn / 2;
1079  constexpr int mid = mm / 2;
1080 
1081  constexpr int stride = Utilities::pow(n_columns, direction);
1082  constexpr int n_blocks1 = stride;
1083  constexpr int n_blocks2 = Utilities::pow(n_rows, (direction>=dim) ? 0 : (dim-direction-1));
1084 
1085  for (int i2=0; i2<n_blocks2; ++i2)
1086  {
1087  for (int i1=0; i1<n_blocks1; ++i1)
1088  {
1089  for (int col=0; col<n_cols; ++col)
1090  {
1091  Number2 val0, val1;
1092  Number in0, in1, res0, res1;
1093  if (contract_over_rows == true)
1094  {
1095  val0 = shape_gradients[col];
1096  val1 = shape_gradients[nn-1-col];
1097  }
1098  else
1099  {
1100  val0 = shape_gradients[col*n_columns];
1101  val1 = shape_gradients[(nn-col-1)*n_columns];
1102  }
1103  if (mid > 0)
1104  {
1105  in0 = in[0];
1106  in1 = in[stride*(mm-1)];
1107  res0 = val0 * in0;
1108  res1 = val1 * in0;
1109  res0 -= val1 * in1;
1110  res1 -= val0 * in1;
1111  for (int ind=1; ind<mid; ++ind)
1112  {
1113  if (contract_over_rows == true)
1114  {
1115  val0 = shape_gradients[ind*n_columns+col];
1116  val1 = shape_gradients[ind*n_columns+nn-1-col];
1117  }
1118  else
1119  {
1120  val0 = shape_gradients[col*n_columns+ind];
1121  val1 = shape_gradients[(nn-col-1)*n_columns+ind];
1122  }
1123  in0 = in[stride*ind];
1124  in1 = in[stride*(mm-1-ind)];
1125  res0 += val0 * in0;
1126  res1 += val1 * in0;
1127  res0 -= val1 * in1;
1128  res1 -= val0 * in1;
1129  }
1130  }
1131  else
1132  res0 = res1 = Number();
1133  if (mm % 2 == 1)
1134  {
1135  if (contract_over_rows == true)
1136  val0 = shape_gradients[mid*n_columns+col];
1137  else
1138  val0 = shape_gradients[col*n_columns+mid];
1139  in1 = val0 * in[stride*mid];
1140  res0 += in1;
1141  res1 -= in1;
1142  }
1143  if (add == false)
1144  {
1145  out[stride*col] = res0;
1146  out[stride*(nn-1-col)] = res1;
1147  }
1148  else
1149  {
1150  out[stride*col] += res0;
1151  out[stride*(nn-1-col)] += res1;
1152  }
1153  }
1154  if ( nn%2 == 1 )
1155  {
1156  Number2 val0;
1157  Number res0;
1158  if (contract_over_rows == true)
1159  val0 = shape_gradients[n_cols];
1160  else
1161  val0 = shape_gradients[n_cols*n_columns];
1162  res0 = val0 * (in[0] - in[stride*(mm-1)]);
1163  for (int ind=1; ind<mid; ++ind)
1164  {
1165  if (contract_over_rows == true)
1166  val0 = shape_gradients[ind*n_columns+n_cols];
1167  else
1168  val0 = shape_gradients[n_cols*n_columns+ind];
1169  Number in1 = val0 * (in[stride*ind] - in[stride*(mm-1-ind)]);
1170  res0 += in1;
1171  }
1172  if (add == false)
1173  out[stride*n_cols] = res0;
1174  else
1175  out[stride*n_cols] += res0;
1176  }
1177 
1178  ++in;
1179  ++out;
1180  }
1181  in += stride*(mm-1);
1182  out += stride*(nn-1);
1183  }
1184  }
1185 
1186 
1187 
1188  // evaluates the given shape data in 1d-3d using the tensor product
1189  // form assuming the symmetries of unit cell shape hessians for
1190  // finite elements in FEEvaluation
1191  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
1192  template <int direction, bool contract_over_rows, bool add>
1193  inline
1194  void
1195  EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
1196  ::hessians (const Number in [],
1197  Number out []) const
1198  {
1199  Assert (shape_hessians != nullptr, ExcNotInitialized());
1200  AssertIndexRange (direction, dim);
1201  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1202  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1203  constexpr int n_cols = nn / 2;
1204  constexpr int mid = mm / 2;
1205 
1206  constexpr int stride = Utilities::pow(n_columns, direction);
1207  constexpr int n_blocks1 = stride;
1208  constexpr int n_blocks2 = Utilities::pow(n_rows, (direction>=dim) ? 0 : (dim-direction-1));
1209 
1210  for (int i2=0; i2<n_blocks2; ++i2)
1211  {
1212  for (int i1=0; i1<n_blocks1; ++i1)
1213  {
1214  for (int col=0; col<n_cols; ++col)
1215  {
1216  Number2 val0, val1;
1217  Number in0, in1, res0, res1;
1218  if (contract_over_rows == true)
1219  {
1220  val0 = shape_hessians[col];
1221  val1 = shape_hessians[nn-1-col];
1222  }
1223  else
1224  {
1225  val0 = shape_hessians[col*n_columns];
1226  val1 = shape_hessians[(col+1)*n_columns-1];
1227  }
1228  if (mid > 0)
1229  {
1230  in0 = in[0];
1231  in1 = in[stride*(mm-1)];
1232  res0 = val0 * in0;
1233  res1 = val1 * in0;
1234  res0 += val1 * in1;
1235  res1 += val0 * in1;
1236  for (int ind=1; ind<mid; ++ind)
1237  {
1238  if (contract_over_rows == true)
1239  {
1240  val0 = shape_hessians[ind*n_columns+col];
1241  val1 = shape_hessians[ind*n_columns+nn-1-col];
1242  }
1243  else
1244  {
1245  val0 = shape_hessians[col*n_columns+ind];
1246  val1 = shape_hessians[(col+1)*n_columns-1-ind];
1247  }
1248  in0 = in[stride*ind];
1249  in1 = in[stride*(mm-1-ind)];
1250  res0 += val0 * in0;
1251  res1 += val1 * in0;
1252  res0 += val1 * in1;
1253  res1 += val0 * in1;
1254  }
1255  }
1256  else
1257  res0 = res1 = Number();
1258  if (mm % 2 == 1)
1259  {
1260  if (contract_over_rows == true)
1261  val0 = shape_hessians[mid*n_columns+col];
1262  else
1263  val0 = shape_hessians[col*n_columns+mid];
1264  in1 = val0 * in[stride*mid];
1265  res0 += in1;
1266  res1 += in1;
1267  }
1268  if (add == false)
1269  {
1270  out[stride*col] = res0;
1271  out[stride*(nn-1-col)] = res1;
1272  }
1273  else
1274  {
1275  out[stride*col] += res0;
1276  out[stride*(nn-1-col)] += res1;
1277  }
1278  }
1279  if ( nn%2 == 1 )
1280  {
1281  Number2 val0;
1282  Number res0;
1283  if (contract_over_rows == true)
1284  val0 = shape_hessians[n_cols];
1285  else
1286  val0 = shape_hessians[n_cols*n_columns];
1287  if (mid > 0)
1288  {
1289  res0 = val0 * (in[0] + in[stride*(mm-1)]);
1290  for (int ind=1; ind<mid; ++ind)
1291  {
1292  if (contract_over_rows == true)
1293  val0 = shape_hessians[ind*n_columns+n_cols];
1294  else
1295  val0 = shape_hessians[n_cols*n_columns+ind];
1296  Number in1 = val0*(in[stride*ind] + in[stride*(mm-1-ind)]);
1297  res0 += in1;
1298  }
1299  }
1300  else
1301  res0 = Number();
1302  if (mm % 2 == 1)
1303  {
1304  if (contract_over_rows == true)
1305  val0 = shape_hessians[mid*n_columns+n_cols];
1306  else
1307  val0 = shape_hessians[n_cols*n_columns+mid];
1308  res0 += val0 * in[stride*mid];
1309  }
1310  if (add == false)
1311  out[stride*n_cols] = res0;
1312  else
1313  out[stride*n_cols] += res0;
1314  }
1315 
1316  ++in;
1317  ++out;
1318  }
1319  in += stride*(mm-1);
1320  out += stride*(nn-1);
1321  }
1322  }
1323 
1324 
1325 
1357  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
1358  struct EvaluatorTensorProduct<evaluate_evenodd,dim,n_rows,n_columns,Number,Number2>
1359  {
1360  static constexpr unsigned int n_rows_of_product = Utilities::pow(n_rows, dim);
1361  static constexpr unsigned int n_columns_of_product = Utilities::pow(n_columns, dim);
1362 
1368  EvaluatorTensorProduct ()
1369  :
1370  shape_values (nullptr),
1371  shape_gradients (nullptr),
1372  shape_hessians (nullptr)
1373  {}
1374 
1379  EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values)
1380  :
1381  shape_values (shape_values.begin()),
1382  shape_gradients (nullptr),
1383  shape_hessians (nullptr)
1384  {
1385  AssertDimension(shape_values.size(), n_rows*((n_columns+1)/2));
1386  }
1387 
1392  EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
1393  const AlignedVector<Number2> &shape_gradients,
1394  const AlignedVector<Number2> &shape_hessians,
1395  const unsigned int dummy1 = 0,
1396  const unsigned int dummy2 = 0)
1397  :
1398  shape_values (shape_values.begin()),
1399  shape_gradients (shape_gradients.begin()),
1400  shape_hessians (shape_hessians.begin())
1401  {
1402  // In this function, we allow for dummy pointers if some of values,
1403  // gradients or hessians should not be computed
1404  if (!shape_values.empty())
1405  AssertDimension(shape_values.size(), n_rows*((n_columns+1)/2));
1406  if (!shape_gradients.empty())
1407  AssertDimension(shape_gradients.size(), n_rows*((n_columns+1)/2));
1408  if (!shape_hessians.empty())
1409  AssertDimension(shape_hessians.size(), n_rows*((n_columns+1)/2));
1410  (void)dummy1;
1411  (void)dummy2;
1412  }
1413 
1414  template <int direction, bool contract_over_rows, bool add>
1415  void
1416  values (const Number in [],
1417  Number out[]) const
1418  {
1419  Assert (shape_values != nullptr, ExcNotInitialized());
1420  apply<direction,contract_over_rows,add,0>(shape_values, in, out);
1421  }
1422 
1423  template <int direction, bool contract_over_rows, bool add>
1424  void
1425  gradients (const Number in [],
1426  Number out[]) const
1427  {
1428  Assert (shape_gradients != nullptr, ExcNotInitialized());
1429  apply<direction,contract_over_rows,add,1>(shape_gradients, in, out);
1430  }
1431 
1432  template <int direction, bool contract_over_rows, bool add>
1433  void
1434  hessians (const Number in [],
1435  Number out[]) const
1436  {
1437  Assert (shape_hessians != nullptr, ExcNotInitialized());
1438  apply<direction,contract_over_rows,add,2>(shape_hessians, in, out);
1439  }
1440 
1441  template <int direction, bool contract_over_rows, bool add>
1442  void
1443  values_one_line (const Number in [],
1444  Number out[]) const
1445  {
1446  Assert (shape_values != nullptr, ExcNotInitialized());
1447  apply<direction,contract_over_rows,add,0,true>(shape_values, in, out);
1448  }
1449 
1450  template <int direction, bool contract_over_rows, bool add>
1451  void
1452  gradients_one_line (const Number in [],
1453  Number out[]) const
1454  {
1455  Assert (shape_gradients != nullptr, ExcNotInitialized());
1456  apply<direction,contract_over_rows,add,1,true>(shape_gradients, in, out);
1457  }
1458 
1459  template <int direction, bool contract_over_rows, bool add>
1460  void
1461  hessians_one_line (const Number in [],
1462  Number out[]) const
1463  {
1464  Assert (shape_hessians != nullptr, ExcNotInitialized());
1465  apply<direction,contract_over_rows,add,2,true>(shape_hessians, in, out);
1466  }
1467 
1496  template <int direction, bool contract_over_rows, bool add, int type,
1497  bool one_line=false>
1498  static void apply (const Number2 *DEAL_II_RESTRICT shape_data,
1499  const Number *in,
1500  Number *out);
1501 
1502  const Number2 *shape_values;
1503  const Number2 *shape_gradients;
1504  const Number2 *shape_hessians;
1505  };
1506 
1507 
1508 
1509  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
1510  template <int direction, bool contract_over_rows, bool add, int type, bool one_line>
1511  inline
1512  void
1513  EvaluatorTensorProduct<evaluate_evenodd,dim,n_rows,n_columns,Number,Number2>
1514  ::apply (const Number2 *DEAL_II_RESTRICT shapes,
1515  const Number *in,
1516  Number *out)
1517  {
1518  static_assert (type < 3, "Only three variants type=0,1,2 implemented");
1519  static_assert (one_line == false || direction==dim-1,
1520  "Single-line evaluation only works for direction=dim-1.");
1521  Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
1522  ExcMessage("In-place operation only supported for "
1523  "n_rows==n_columns or single-line interpolation"));
1524 
1525  // We cannot statically assert that direction is less than dim, so must do
1526  // an additional dynamic check
1527  AssertIndexRange (direction, dim);
1528 
1529  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1530  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1531  constexpr int n_cols = nn / 2;
1532  constexpr int mid = mm / 2;
1533 
1534  constexpr int stride = Utilities::pow(n_columns, direction);
1535  constexpr int n_blocks1 = one_line ? 1 : stride;
1536  constexpr int n_blocks2 = Utilities::pow(n_rows, (direction>=dim) ? 0 : (dim-direction-1));
1537 
1538  constexpr int offset = (n_columns+1)/2;
1539 
1540  // this code may look very inefficient at first sight due to the many
1541  // different cases with if's at the innermost loop part, but all of the
1542  // conditionals can be evaluated at compile time because they are
1543  // templates, so the compiler should optimize everything away
1544  for (int i2=0; i2<n_blocks2; ++i2)
1545  {
1546  for (int i1=0; i1<n_blocks1; ++i1)
1547  {
1548  Number xp[mid>0?mid:1], xm[mid>0?mid:1];
1549  for (int i=0; i<mid; ++i)
1550  {
1551  if (contract_over_rows == true && type == 1)
1552  {
1553  xp[i] = in[stride*i] - in[stride*(mm-1-i)];
1554  xm[i] = in[stride*i] + in[stride*(mm-1-i)];
1555  }
1556  else
1557  {
1558  xp[i] = in[stride*i] + in[stride*(mm-1-i)];
1559  xm[i] = in[stride*i] - in[stride*(mm-1-i)];
1560  }
1561  }
1562  Number xmid = in[stride*mid];
1563  for (int col=0; col<n_cols; ++col)
1564  {
1565  Number r0, r1;
1566  if (mid > 0)
1567  {
1568  if (contract_over_rows == true)
1569  {
1570  r0 = shapes[col] * xp[0];
1571  r1 = shapes[(n_rows-1)*offset + col] * xm[0];
1572  }
1573  else
1574  {
1575  r0 = shapes[col*offset] * xp[0];
1576  r1 = shapes[(n_rows-1-col)*offset] * xm[0];
1577  }
1578  for (int ind=1; ind<mid; ++ind)
1579  {
1580  if (contract_over_rows == true)
1581  {
1582  r0 += shapes[ind*offset+col] * xp[ind];
1583  r1 += shapes[(n_rows-1-ind)*offset+col] * xm[ind];
1584  }
1585  else
1586  {
1587  r0 += shapes[col*offset+ind] * xp[ind];
1588  r1 += shapes[(n_rows-1-col)*offset+ind] * xm[ind];
1589  }
1590  }
1591  }
1592  else
1593  r0 = r1 = Number();
1594  if (mm % 2 == 1 && contract_over_rows == true)
1595  {
1596  if (type == 1)
1597  r1 += shapes[mid*offset+col] * xmid;
1598  else
1599  r0 += shapes[mid*offset+col] * xmid;
1600  }
1601  else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
1602  r0 += shapes[col*offset+mid] * xmid;
1603 
1604  if (add == false)
1605  {
1606  out[stride*col] = r0 + r1;
1607  if (type == 1 && contract_over_rows == false)
1608  out[stride*(nn-1-col)] = r1 - r0;
1609  else
1610  out[stride*(nn-1-col)] = r0 - r1;
1611  }
1612  else
1613  {
1614  out[stride*col] += r0 + r1;
1615  if (type == 1 && contract_over_rows == false)
1616  out[stride*(nn-1-col)] += r1 - r0;
1617  else
1618  out[stride*(nn-1-col)] += r0 - r1;
1619  }
1620  }
1621  if ( type == 0 && contract_over_rows == true && nn%2==1 && mm%2==1 )
1622  {
1623  if (add==false)
1624  out[stride*n_cols] = shapes[mid*offset+n_cols] * xmid;
1625  else
1626  out[stride*n_cols] += shapes[mid*offset+n_cols] * xmid;
1627  }
1628  else if (contract_over_rows == true && nn%2==1)
1629  {
1630  Number r0;
1631  if (mid > 0)
1632  {
1633  r0 = shapes[n_cols] * xp[0];
1634  for (int ind=1; ind<mid; ++ind)
1635  r0 += shapes[ind*offset+n_cols] * xp[ind];
1636  }
1637  else
1638  r0 = Number();
1639  if (type != 1 && mm % 2 == 1)
1640  r0 += shapes[mid*offset+n_cols] * xmid;
1641 
1642  if (add == false)
1643  out[stride*n_cols] = r0;
1644  else
1645  out[stride*n_cols] += r0;
1646  }
1647  else if (contract_over_rows == false && nn%2 == 1)
1648  {
1649  Number r0;
1650  if (mid > 0)
1651  {
1652  if (type == 1)
1653  {
1654  r0 = shapes[n_cols*offset] * xm[0];
1655  for (int ind=1; ind<mid; ++ind)
1656  r0 += shapes[n_cols*offset+ind] * xm[ind];
1657  }
1658  else
1659  {
1660  r0 = shapes[n_cols*offset] * xp[0];
1661  for (int ind=1; ind<mid; ++ind)
1662  r0 += shapes[n_cols*offset+ind] * xp[ind];
1663  }
1664  }
1665  else
1666  r0 = Number();
1667 
1668  if ((type == 0 || type == 2) && mm % 2 == 1)
1669  r0 += shapes[n_cols*offset+mid] * xmid;
1670 
1671  if (add == false)
1672  out[stride*n_cols] = r0;
1673  else
1674  out[stride*n_cols] += r0;
1675  }
1676  if (one_line == false)
1677  {
1678  in += 1;
1679  out += 1;
1680  }
1681  }
1682  if (one_line == false)
1683  {
1684  in += stride * (mm-1);
1685  out += stride * (nn-1);
1686  }
1687  }
1688  }
1689 
1690 
1691 
1720  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
1721  struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,dim,n_rows,n_columns,Number,Number2>
1722  {
1723  static constexpr unsigned int n_rows_of_product = Utilities::pow(n_rows, dim);
1724  static constexpr unsigned int n_columns_of_product = Utilities::pow(n_columns, dim);
1725 
1731  EvaluatorTensorProduct ()
1732  :
1733  shape_values (nullptr),
1734  shape_gradients (nullptr),
1735  shape_hessians (nullptr)
1736  {}
1737 
1742  EvaluatorTensorProduct (const AlignedVector<Number> &shape_values)
1743  :
1744  shape_values (shape_values.begin()),
1745  shape_gradients (nullptr),
1746  shape_hessians (nullptr)
1747  {}
1748 
1753  EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
1754  const AlignedVector<Number2> &shape_gradients,
1755  const AlignedVector<Number2> &shape_hessians,
1756  const unsigned int dummy1 = 0,
1757  const unsigned int dummy2 = 0)
1758  :
1759  shape_values (shape_values.begin()),
1760  shape_gradients (shape_gradients.begin()),
1761  shape_hessians (shape_hessians.begin())
1762  {
1763  (void)dummy1;
1764  (void)dummy2;
1765  }
1766 
1767  template <int direction, bool contract_over_rows, bool add>
1768  void
1769  values (const Number in [],
1770  Number out[]) const
1771  {
1772  Assert (shape_values != nullptr, ExcNotInitialized());
1773  apply<direction,contract_over_rows,add,0>(shape_values, in, out);
1774  }
1775 
1776  template <int direction, bool contract_over_rows, bool add>
1777  void
1778  gradients (const Number in [],
1779  Number out[]) const
1780  {
1781  Assert (shape_gradients != nullptr, ExcNotInitialized());
1782  apply<direction,contract_over_rows,add,1>(shape_gradients, in, out);
1783  }
1784 
1785  template <int direction, bool contract_over_rows, bool add>
1786  void
1787  hessians (const Number in [],
1788  Number out[]) const
1789  {
1790  Assert (shape_hessians != nullptr, ExcNotInitialized());
1791  apply<direction,contract_over_rows,add,0>(shape_hessians, in, out);
1792  }
1793 
1794  template <int direction, bool contract_over_rows, bool add>
1795  void
1796  values_one_line (const Number in [],
1797  Number out[]) const
1798  {
1799  Assert (shape_values != nullptr, ExcNotInitialized());
1800  apply<direction,contract_over_rows,add,0,true>(shape_values, in, out);
1801  }
1802 
1803  template <int direction, bool contract_over_rows, bool add>
1804  void
1805  gradients_one_line (const Number in [],
1806  Number out[]) const
1807  {
1808  Assert (shape_gradients != nullptr, ExcNotInitialized());
1809  apply<direction,contract_over_rows,add,1,true>(shape_gradients, in, out);
1810  }
1811 
1812  template <int direction, bool contract_over_rows, bool add>
1813  void
1814  hessians_one_line (const Number in [],
1815  Number out[]) const
1816  {
1817  Assert (shape_hessians != nullptr, ExcNotInitialized());
1818  apply<direction,contract_over_rows,add,0,true>(shape_hessians, in, out);
1819  }
1820 
1848  template <int direction, bool contract_over_rows, bool add, int type,
1849  bool one_line=false>
1850  static void apply (const Number2 *DEAL_II_RESTRICT shape_data,
1851  const Number *in,
1852  Number *out);
1853 
1854  const Number2 *shape_values;
1855  const Number2 *shape_gradients;
1856  const Number2 *shape_hessians;
1857  };
1858 
1859 
1860 
1861  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
1862  template <int direction, bool contract_over_rows, bool add, int type, bool one_line>
1863  inline
1864  void
1865  EvaluatorTensorProduct<evaluate_symmetric_hierarchical,dim,n_rows,n_columns,Number,Number2>
1866  ::apply (const Number2 *DEAL_II_RESTRICT shapes,
1867  const Number *in,
1868  Number *out)
1869  {
1870  static_assert (one_line == false || direction==dim-1,
1871  "Single-line evaluation only works for direction=dim-1.");
1872  static_assert (type == 0 || type == 1,
1873  "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
1874  Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
1875  ExcMessage("In-place operation only supported for "
1876  "n_rows==n_columns or single-line interpolation"));
1877 
1878  // We cannot statically assert that direction is less than dim, so must do
1879  // an additional dynamic check
1880  AssertIndexRange (direction, dim);
1881 
1882  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1883  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1884  constexpr int n_cols = nn / 2;
1885  constexpr int mid = mm / 2;
1886 
1887  constexpr int stride = Utilities::pow(n_columns, direction);
1888  constexpr int n_blocks1 = one_line ? 1 : stride;
1889  constexpr int n_blocks2 = Utilities::pow(n_rows, (direction>=dim) ? 0 : (dim-direction-1));
1890 
1891  // this code may look very inefficient at first sight due to the many
1892  // different cases with if's at the innermost loop part, but all of the
1893  // conditionals can be evaluated at compile time because they are
1894  // templates, so the compiler should optimize everything away
1895  for (int i2=0; i2<n_blocks2; ++i2)
1896  {
1897  for (int i1=0; i1<n_blocks1; ++i1)
1898  {
1899  if (contract_over_rows)
1900  {
1901  Number x[mm];
1902  for (unsigned int i=0; i<mm; ++i)
1903  x[i] = in[stride*i];
1904  for (unsigned int col=0; col<n_cols; ++col)
1905  {
1906  Number r0, r1;
1907  if (mid > 0)
1908  {
1909  r0 = shapes[col] * x[0];
1910  r1 = shapes[col+n_columns] * x[1];
1911  for (unsigned int ind=1; ind<mid; ++ind)
1912  {
1913  r0 += shapes[col+2*ind*n_columns] * x[2*ind];
1914  r1 += shapes[col+(2*ind+1)*n_columns] * x[2*ind+1];
1915  }
1916  }
1917  else
1918  r0 = r1 = Number();
1919  if (mm%2 == 1)
1920  r0 += shapes[col+(mm-1)*n_columns] * x[mm-1];
1921  if (add == false)
1922  {
1923  out[stride*col] = r0 + r1;
1924  if (type == 1)
1925  out[stride*(nn-1-col)] = r1 - r0;
1926  else
1927  out[stride*(nn-1-col)] = r0 - r1;
1928  }
1929  else
1930  {
1931  out[stride*col] += r0 + r1;
1932  if (type == 1)
1933  out[stride*(nn-1-col)] += r1 - r0;
1934  else
1935  out[stride*(nn-1-col)] += r0 - r1;
1936  }
1937  }
1938  if (nn%2 == 1)
1939  {
1940  Number r0;
1941  const unsigned int shift = type==1 ? 1 : 0;
1942  if (mid>0)
1943  {
1944  r0 = shapes[n_cols + shift*n_columns] * x[shift];
1945  for (unsigned int ind=1; ind<mid; ++ind)
1946  r0 += shapes[n_cols + (2*ind+shift)*n_columns] * x[2*ind+shift];
1947  }
1948  else
1949  r0 = 0;
1950  if (type != 1 && mm%2 == 1)
1951  r0 += shapes[n_cols + (mm-1)*n_columns] * x[mm-1];
1952  if (add == false)
1953  out[stride*n_cols] = r0;
1954  else
1955  out[stride*n_cols] += r0;
1956  }
1957  }
1958  else
1959  {
1960  Number xp[mid+1], xm[mid>0?mid:1];
1961  for (int i=0; i<mid; ++i)
1962  if (type == 0)
1963  {
1964  xp[i] = in[stride*i] + in[stride*(mm-1-i)];
1965  xm[i] = in[stride*i] - in[stride*(mm-1-i)];
1966  }
1967  else
1968  {
1969  xp[i] = in[stride*i] - in[stride*(mm-1-i)];
1970  xm[i] = in[stride*i] + in[stride*(mm-1-i)];
1971  }
1972  if (mm%2 == 1)
1973  xp[mid] = in[stride*mid];
1974  for (unsigned int col=0; col<n_cols; ++col)
1975  {
1976  Number r0, r1;
1977  if (mid > 0)
1978  {
1979  r0 = shapes[2*col*n_columns] * xp[0];
1980  r1 = shapes[(2*col+1)*n_columns] * xm[0];
1981  for (unsigned int ind=1; ind<mid; ++ind)
1982  {
1983  r0 += shapes[2*col*n_columns+ind] * xp[ind];
1984  r1 += shapes[(2*col+1)*n_columns+ind] * xm[ind];
1985  }
1986  }
1987  else
1988  r0 = r1 = Number();
1989  if (mm%2 == 1)
1990  {
1991  if (type == 1)
1992  r1 += shapes[(2*col+1)*n_columns+mid] * xp[mid];
1993  else
1994  r0 += shapes[2*col*n_columns+mid] * xp[mid];
1995  }
1996  if (add == false)
1997  {
1998  out[stride*(2*col)] = r0;
1999  out[stride*(2*col+1)] = r1;
2000  }
2001  else
2002  {
2003  out[stride*(2*col)] += r0;
2004  out[stride*(2*col+1)] += r1;
2005  }
2006  }
2007  if (nn%2 == 1)
2008  {
2009  Number r0;
2010  if (mid > 0)
2011  {
2012  r0 = shapes[(nn-1)*n_columns] * xp[0];
2013  for (unsigned int ind=1; ind<mid; ++ind)
2014  r0 += shapes[(nn-1)*n_columns+ind] * xp[ind];
2015  }
2016  else
2017  r0 = Number();
2018  if (mm%2 == 1 && type == 0)
2019  r0 += shapes[(nn-1)*n_columns+mid] * xp[mid];
2020  if (add == false)
2021  out[stride*(nn-1)] = r0;
2022  else
2023  out[stride*(nn-1)] += r0;
2024  }
2025  }
2026  if (one_line == false)
2027  {
2028  in += 1;
2029  out += 1;
2030  }
2031  }
2032  if (one_line == false)
2033  {
2034  in += stride * (mm-1);
2035  out += stride * (nn-1);
2036  }
2037  }
2038  }
2039 
2040 } // end of namespace internal
2041 
2042 
2043 DEAL_II_NAMESPACE_CLOSE
2044 
2045 #endif
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
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::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:1392
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1284
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition: tensor_product_kernels.h:1514
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:863
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number > &shape_values)
Definition: tensor_product_kernels.h:1742
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:1753
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1142
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:1731
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:123
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
Definition: tensor_product_kernels.h:509
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:133
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values)
Definition: tensor_product_kernels.h:1379
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:825
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:497
AlignedVector::size
size_type size() const
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:1368
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition: tensor_product_kernels.h:286
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
AlignedVector::empty
bool empty() const
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition: tensor_product_kernels.h:1866
Generated by   doxygen 1.8.14