Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor_product_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_tensor_product_kernels_h
18 #define dealii_matrix_free_tensor_product_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/ndarray.h>
25 #include <deal.II/base/utilities.h>
26 
28 
29 
31 
32 
33 
34 namespace internal
35 {
41  {
70  };
71 
72 
73 
77  enum class EvaluatorQuantity
78  {
82  value,
86  gradient,
90  hessian
91  };
92 
93 
94 
110  template <EvaluatorVariant variant,
111  EvaluatorQuantity quantity,
112  int n_rows,
113  int n_columns,
114  int stride_in,
115  int stride_out,
116  bool transpose_matrix,
117  bool add,
118  typename Number,
119  typename Number2>
120  std::enable_if_t<(variant == evaluate_general), void>
122  const Number *in,
123  Number *out)
124  {
125  // We can only statically assert that one argument is non-zero because
126  // face evaluation might instantiate some functions, so we need to use the
127  // run-time assert to verify that we do not end up involuntarily.
128  static_assert(n_rows > 0 || n_columns > 0,
129  "Specialization only for n_rows, n_columns > 0");
130  Assert(n_rows > 0 && n_columns > 0,
131  ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
132  std::to_string(n_rows) + ", " +
133  std::to_string(n_columns) + " was passed!"));
134  static_assert(quantity == EvaluatorQuantity::value,
135  "This function should only use EvaluatorQuantity::value");
136 
137  constexpr int mm = transpose_matrix ? n_rows : n_columns,
138  nn = transpose_matrix ? n_columns : n_rows;
139 
140  std::array<Number, mm> x;
141  for (int i = 0; i < mm; ++i)
142  x[i] = in[stride_in * i];
143  for (int col = 0; col < nn; ++col)
144  {
145  Number res0;
146  if (transpose_matrix == true)
147  {
148  res0 = matrix[col] * x[0];
149  for (int i = 1; i < mm; ++i)
150  res0 += matrix[i * n_columns + col] * x[i];
151  }
152  else
153  {
154  res0 = matrix[col * n_columns] * x[0];
155  for (int i = 1; i < mm; ++i)
156  res0 += matrix[col * n_columns + i] * x[i];
157  }
158  if (add)
159  out[stride_out * col] += res0;
160  else
161  out[stride_out * col] = res0;
162  }
163  }
164 
165 
166 
171  template <EvaluatorVariant variant,
172  EvaluatorQuantity quantity,
173  bool transpose_matrix,
174  bool add,
175  typename Number,
176  typename Number2>
177  std::enable_if_t<(variant == evaluate_general), void>
179  const Number *in,
180  Number *out,
181  const int n_rows,
182  const int n_columns,
183  const int stride_in,
184  const int stride_out)
185  {
186  const int mm = transpose_matrix ? n_rows : n_columns,
187  nn = transpose_matrix ? n_columns : n_rows;
188  Assert(n_rows <= 128, ExcNotImplemented());
189  Assert(n_rows > 0 && n_columns > 0,
190  ExcInternalError("Empty evaluation task!"));
191  Assert(n_rows > 0 && n_columns > 0,
192  ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
193  std::to_string(n_rows) + ", " +
194  std::to_string(n_columns) + " was passed!"));
195 
196  static_assert(quantity == EvaluatorQuantity::value,
197  "This function should only use EvaluatorQuantity::value");
198 
199  // specialization for n_rows = 2 that manually unrolls the innermost loop
200  // to make the operation perform better (not completely as good as the
201  // templated one, but much better than the generic version down below,
202  // because the loop over col can be more effectively unrolled by the
203  // compiler)
204  if (transpose_matrix && n_rows == 2)
205  {
206  const Number2 *matrix_1 = matrix + n_columns;
207  const Number x0 = in[0], x1 = in[stride_in];
208  for (int col = 0; col < nn; ++col)
209  {
210  const Number result = matrix[col] * x0 + matrix_1[col] * x1;
211  if (add)
212  out[stride_out * col] += result;
213  else
214  out[stride_out * col] = result;
215  }
216  }
217  else if (transpose_matrix && n_rows == 3)
218  {
219  const Number2 *matrix_1 = matrix + n_columns;
220  const Number2 *matrix_2 = matrix_1 + n_columns;
221  const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
222  for (int col = 0; col < nn; ++col)
223  {
224  const Number result =
225  matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
226  if (add)
227  out[stride_out * col] += result;
228  else
229  out[stride_out * col] = result;
230  }
231  }
232  else
233  {
234  std::array<Number, 129> x;
235  for (int i = 0; i < mm; ++i)
236  x[i] = in[stride_in * i];
237 
238  Number res0;
239  for (int col = 0; col < nn; ++col)
240  {
241  if (transpose_matrix == true)
242  {
243  res0 = matrix[col] * x[0];
244  for (int i = 1; i < mm; ++i)
245  res0 += matrix[i * n_columns + col] * x[i];
246  }
247  else
248  {
249  res0 = matrix[col * n_columns] * x[0];
250  for (int i = 1; i < mm; ++i)
251  res0 += matrix[col * n_columns + i] * x[i];
252  }
253  if (add)
254  out[stride_out * col] += res0;
255  else
256  out[stride_out * col] = res0;
257  }
258  }
259  }
260 
261 
262 
269  template <EvaluatorVariant variant,
270  EvaluatorQuantity quantity,
271  int n_rows,
272  int n_columns,
273  int stride_in,
274  int stride_out,
275  bool transpose_matrix,
276  bool add,
277  typename Number,
278  typename Number2>
279  std::enable_if_t<(variant == evaluate_symmetric), void>
281  const Number *in,
282  Number *out)
283  {
284  // We can only statically assert that one argument is non-zero because
285  // face evaluation might instantiate some functions, so we need to use the
286  // run-time assert to verify that we do not end up involuntarily.
287  static_assert(n_rows > 0 || n_columns > 0,
288  "Specialization only for n_rows, n_columns > 0");
289  Assert(n_rows > 0 && n_columns > 0,
290  ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
291  std::to_string(n_rows) + ", " +
292  std::to_string(n_columns) + " was passed!"));
293 
294  constexpr int mm = transpose_matrix ? n_rows : n_columns,
295  nn = transpose_matrix ? n_columns : n_rows;
296  constexpr int n_cols = nn / 2;
297  constexpr int mid = mm / 2;
298 
299  std::array<Number, mm> x;
300  for (int i = 0; i < mm; ++i)
301  x[i] = in[stride_in * i];
302 
303  if (quantity == EvaluatorQuantity::value)
304  {
305  // In this case, the 1d shape values read (sorted lexicographically,
306  // rows run over 1d dofs, columns over quadrature points):
307  // Q2 --> [ 0.687 0 -0.087 ]
308  // [ 0.4 1 0.4 ]
309  // [-0.087 0 0.687 ]
310  // Q3 --> [ 0.66 0.003 0.002 0.049 ]
311  // [ 0.521 1.005 -0.01 -0.230 ]
312  // [-0.230 -0.01 1.005 0.521 ]
313  // [ 0.049 0.002 0.003 0.66 ]
314  // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
315  // [ 0.608 1.059 0 0.039 0.176 ]
316  // [-0.409 -0.113 1 -0.113 -0.409 ]
317  // [ 0.176 0.039 0 1.059 0.608 ]
318  // [-0.032 -0.007 0 0.022 0.658 ]
319  //
320  // In these matrices, we want to use avoid computations involving
321  // zeros and ones and use the symmetry in entries starting from (1,1)
322  // forward and (N,N) backward, respectively to reduce the number of
323  // read operations.
324  for (int col = 0; col < n_cols; ++col)
325  {
326  Number2 val0, val1;
327  Number res0, res1;
328  if (transpose_matrix == true)
329  {
330  val0 = matrix[col];
331  val1 = matrix[nn - 1 - col];
332  }
333  else
334  {
335  val0 = matrix[col * n_columns];
336  val1 = matrix[(col + 1) * n_columns - 1];
337  }
338  if (mid > 0)
339  {
340  res0 = val0 * x[0];
341  res1 = val1 * x[0];
342  res0 += val1 * x[mm - 1];
343  res1 += val0 * x[mm - 1];
344  for (int ind = 1; ind < mid; ++ind)
345  {
346  if (transpose_matrix == true)
347  {
348  val0 = matrix[ind * n_columns + col];
349  val1 = matrix[ind * n_columns + nn - 1 - col];
350  }
351  else
352  {
353  val0 = matrix[col * n_columns + ind];
354  val1 = matrix[(col + 1) * n_columns - 1 - ind];
355  }
356  res0 += val0 * x[ind];
357  res1 += val1 * x[ind];
358  res0 += val1 * x[mm - 1 - ind];
359  res1 += val0 * x[mm - 1 - ind];
360  }
361  }
362  else
363  res0 = res1 = Number();
364  if (transpose_matrix == true)
365  {
366  if (mm % 2 == 1)
367  {
368  const Number tmp = matrix[mid * n_columns + col] * x[mid];
369  res0 += tmp;
370  res1 += tmp;
371  }
372  }
373  else
374  {
375  if (mm % 2 == 1 && nn % 2 == 0)
376  {
377  const Number tmp = matrix[col * n_columns + mid] * x[mid];
378  res0 += tmp;
379  res1 += tmp;
380  }
381  }
382  if (add)
383  {
384  out[stride_out * col] += res0;
385  out[stride_out * (nn - 1 - col)] += res1;
386  }
387  else
388  {
389  out[stride_out * col] = res0;
390  out[stride_out * (nn - 1 - col)] = res1;
391  }
392  }
393  if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1)
394  {
395  if (add)
396  out[stride_out * n_cols] += x[mid];
397  else
398  out[stride_out * n_cols] = x[mid];
399  }
400  else if (transpose_matrix == true && nn % 2 == 1)
401  {
402  Number res0;
403  if (mid > 0)
404  {
405  res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
406  for (int ind = 1; ind < mid; ++ind)
407  {
408  const Number2 val0 = matrix[ind * n_columns + n_cols];
409  res0 += val0 * (x[ind] + in[mm - 1 - ind]);
410  }
411  }
412  else
413  res0 = Number();
414  if (add)
415  out[stride_out * n_cols] += res0;
416  else
417  out[stride_out * n_cols] = res0;
418  }
419  else if (transpose_matrix == false && nn % 2 == 1)
420  {
421  Number res0;
422  if (mid > 0)
423  {
424  res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
425  for (int ind = 1; ind < mid; ++ind)
426  {
427  const Number2 val0 = matrix[n_cols * n_columns + ind];
428  res0 += val0 * (x[ind] + x[mm - 1 - ind]);
429  ;
430  }
431  if (mm % 2)
432  res0 += x[mid];
433  }
434  else
435  res0 = in[0];
436  if (add)
437  out[stride_out * n_cols] += res0;
438  else
439  out[stride_out * n_cols] = res0;
440  }
441  }
442  else if (quantity == EvaluatorQuantity::gradient)
443  {
444  // For the specialized loop used for gradient computations we again
445  // exploit symmetries according to the following entries (sorted
446  // lexicographically, rows run over 1d dofs, columns over quadrature
447  // points):
448  // Q2 --> [-2.549 -1 0.549 ]
449  // [ 3.098 0 -3.098 ]
450  // [-0.549 1 2.549 ]
451  // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
452  // [ 6.07 -1.44 -2.97 2.196 ]
453  // [-2.196 2.97 1.44 -6.07 ]
454  // [ 0.44 -0.5 1.03 4.315 ]
455  // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
456  // [10.111 -2.76 -2.667 2.066 -2.306 ]
457  // [-5.688 5.773 0 -5.773 5.688 ]
458  // [ 2.306 -2.066 2.667 2.76 -10.111 ]
459  // [-0.413 0.353 -0.333 -0.353 0.413 ]
460  for (int col = 0; col < n_cols; ++col)
461  {
462  Number2 val0, val1;
463  Number res0, res1;
464  if (transpose_matrix == true)
465  {
466  val0 = matrix[col];
467  val1 = matrix[nn - 1 - col];
468  }
469  else
470  {
471  val0 = matrix[col * n_columns];
472  val1 = matrix[(nn - col - 1) * n_columns];
473  }
474  if (mid > 0)
475  {
476  res0 = val0 * x[0];
477  res1 = val1 * x[0];
478  res0 -= val1 * x[mm - 1];
479  res1 -= val0 * x[mm - 1];
480  for (int ind = 1; ind < mid; ++ind)
481  {
482  if (transpose_matrix == true)
483  {
484  val0 = matrix[ind * n_columns + col];
485  val1 = matrix[ind * n_columns + nn - 1 - col];
486  }
487  else
488  {
489  val0 = matrix[col * n_columns + ind];
490  val1 = matrix[(nn - col - 1) * n_columns + ind];
491  }
492  res0 += val0 * x[ind];
493  res1 += val1 * x[ind];
494  res0 -= val1 * x[mm - 1 - ind];
495  res1 -= val0 * x[mm - 1 - ind];
496  }
497  }
498  else
499  res0 = res1 = Number();
500  if (mm % 2 == 1)
501  {
502  if (transpose_matrix == true)
503  val0 = matrix[mid * n_columns + col];
504  else
505  val0 = matrix[col * n_columns + mid];
506  const Number tmp = val0 * x[mid];
507  res0 += tmp;
508  res1 -= tmp;
509  }
510  if (add)
511  {
512  out[stride_out * col] += res0;
513  out[stride_out * (nn - 1 - col)] += res1;
514  }
515  else
516  {
517  out[stride_out * col] = res0;
518  out[stride_out * (nn - 1 - col)] = res1;
519  }
520  }
521  if (nn % 2 == 1)
522  {
523  Number2 val0;
524  Number res0;
525  if (transpose_matrix == true)
526  val0 = matrix[n_cols];
527  else
528  val0 = matrix[n_cols * n_columns];
529  res0 = val0 * (x[0] - x[mm - 1]);
530  for (int ind = 1; ind < mid; ++ind)
531  {
532  if (transpose_matrix == true)
533  val0 = matrix[ind * n_columns + n_cols];
534  else
535  val0 = matrix[n_cols * n_columns + ind];
536  Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
537  res0 += in1;
538  }
539  if (add)
540  out[stride_out * n_cols] += res0;
541  else
542  out[stride_out * n_cols] = res0;
543  }
544  }
545  else
546  {
547  // Hessians are almost the same as values, apart from some missing '1'
548  // entries
549  for (int col = 0; col < n_cols; ++col)
550  {
551  Number2 val0, val1;
552  Number res0, res1;
553  if (transpose_matrix == true)
554  {
555  val0 = matrix[col];
556  val1 = matrix[nn - 1 - col];
557  }
558  else
559  {
560  val0 = matrix[col * n_columns];
561  val1 = matrix[(col + 1) * n_columns - 1];
562  }
563  if (mid > 0)
564  {
565  res0 = val0 * x[0];
566  res1 = val1 * x[0];
567  res0 += val1 * x[mm - 1];
568  res1 += val0 * x[mm - 1];
569  for (int ind = 1; ind < mid; ++ind)
570  {
571  if (transpose_matrix == true)
572  {
573  val0 = matrix[ind * n_columns + col];
574  val1 = matrix[ind * n_columns + nn - 1 - col];
575  }
576  else
577  {
578  val0 = matrix[col * n_columns + ind];
579  val1 = matrix[(col + 1) * n_columns - 1 - ind];
580  }
581  res0 += val0 * x[ind];
582  res1 += val1 * x[ind];
583  res0 += val1 * x[mm - 1 - ind];
584  res1 += val0 * x[mm - 1 - ind];
585  }
586  }
587  else
588  res0 = res1 = Number();
589  if (mm % 2 == 1)
590  {
591  if (transpose_matrix == true)
592  val0 = matrix[mid * n_columns + col];
593  else
594  val0 = matrix[col * n_columns + mid];
595  const Number tmp = val0 * x[mid];
596  res0 += tmp;
597  res1 += tmp;
598  }
599  if (add)
600  {
601  out[stride_out * col] += res0;
602  out[stride_out * (nn - 1 - col)] += res1;
603  }
604  else
605  {
606  out[stride_out * col] = res0;
607  out[stride_out * (nn - 1 - col)] = res1;
608  }
609  }
610  if (nn % 2 == 1)
611  {
612  Number2 val0;
613  Number res0;
614  if (transpose_matrix == true)
615  val0 = matrix[n_cols];
616  else
617  val0 = matrix[n_cols * n_columns];
618  if (mid > 0)
619  {
620  res0 = val0 * (x[0] + x[mm - 1]);
621  for (int ind = 1; ind < mid; ++ind)
622  {
623  if (transpose_matrix == true)
624  val0 = matrix[ind * n_columns + n_cols];
625  else
626  val0 = matrix[n_cols * n_columns + ind];
627  Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
628  res0 += in1;
629  }
630  }
631  else
632  res0 = Number();
633  if (mm % 2 == 1)
634  {
635  if (transpose_matrix == true)
636  val0 = matrix[mid * n_columns + n_cols];
637  else
638  val0 = matrix[n_cols * n_columns + mid];
639  res0 += val0 * x[mid];
640  }
641  if (add)
642  out[stride_out * n_cols] += res0;
643  else
644  out[stride_out * n_cols] = res0;
645  }
646  }
647  }
648 
649 
650 
669  template <EvaluatorVariant variant,
670  EvaluatorQuantity quantity,
671  int n_rows_static,
672  int n_columns_static,
673  int stride_in_static,
674  int stride_out_static,
675  bool transpose_matrix,
676  bool add,
677  typename Number,
678  typename Number2>
679 #ifndef DEBUG
680  inline DEAL_II_ALWAYS_INLINE
681 #endif
682  std::enable_if_t<(variant == evaluate_evenodd), void>
684  const Number *in,
685  Number *out,
686  int n_rows_runtime = 0,
687  int n_columns_runtime = 0,
688  int stride_in_runtime = 0,
689  int stride_out_runtime = 0)
690  {
691  static_assert(n_rows_static >= 0 && n_columns_static >= 0,
692  "Negative loop ranges are not allowed!");
693 
694  const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
695  const int n_columns =
696  n_rows_static == 0 ? n_columns_runtime : n_columns_static;
697  const int stride_in =
698  n_rows_static == 0 ? stride_in_runtime : stride_in_static;
699  const int stride_out =
700  n_rows_static == 0 ? stride_out_runtime : stride_out_static;
701 
702  Assert(n_rows > 0 && n_columns > 0,
703  ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
704  std::to_string(n_rows) + ", " +
705  std::to_string(n_columns) + " was passed!"));
706 
707  const int mm = transpose_matrix ? n_rows : n_columns,
708  nn = transpose_matrix ? n_columns : n_rows;
709  const int n_half = nn / 2;
710  const int m_half = mm / 2;
711 
712  constexpr int array_length =
713  (n_rows_static == 0) ?
714  16 // for non-templated execution
715  :
716  (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
717  const int offset = (n_columns + 1) / 2;
718 
719  Assert(m_half <= array_length, ExcNotImplemented());
720 
721  std::array<Number, array_length> xp, xm;
722  for (int i = 0; i < m_half; ++i)
723  {
724  if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient)
725  {
726  xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
727  xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
728  }
729  else
730  {
731  xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
732  xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
733  }
734  }
735  Number xmid = in[stride_in * m_half];
736  for (int col = 0; col < n_half; ++col)
737  {
738  Number r0, r1;
739  if (m_half > 0)
740  {
741  if (transpose_matrix == true)
742  {
743  r0 = matrix[col] * xp[0];
744  r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
745  }
746  else
747  {
748  r0 = matrix[col * offset] * xp[0];
749  r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
750  }
751  for (int ind = 1; ind < m_half; ++ind)
752  {
753  if (transpose_matrix == true)
754  {
755  r0 += matrix[ind * offset + col] * xp[ind];
756  r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
757  }
758  else
759  {
760  r0 += matrix[col * offset + ind] * xp[ind];
761  r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
762  }
763  }
764  }
765  else
766  r0 = r1 = Number();
767  if (mm % 2 == 1 && transpose_matrix == true)
768  {
769  if (quantity == EvaluatorQuantity::gradient)
770  r1 += matrix[m_half * offset + col] * xmid;
771  else
772  r0 += matrix[m_half * offset + col] * xmid;
773  }
774  else if (mm % 2 == 1 &&
775  (nn % 2 == 0 || quantity != EvaluatorQuantity::value ||
776  mm == 3))
777  r0 += matrix[col * offset + m_half] * xmid;
778 
779  if (add)
780  {
781  out[stride_out * col] += r0 + r1;
782  if (quantity == EvaluatorQuantity::gradient &&
783  transpose_matrix == false)
784  out[stride_out * (nn - 1 - col)] += r1 - r0;
785  else
786  out[stride_out * (nn - 1 - col)] += r0 - r1;
787  }
788  else
789  {
790  out[stride_out * col] = r0 + r1;
791  if (quantity == EvaluatorQuantity::gradient &&
792  transpose_matrix == false)
793  out[stride_out * (nn - 1 - col)] = r1 - r0;
794  else
795  out[stride_out * (nn - 1 - col)] = r0 - r1;
796  }
797  }
798  if (quantity == EvaluatorQuantity::value && transpose_matrix == true &&
799  nn % 2 == 1 && mm % 2 == 1 && mm > 3)
800  {
801  if (add)
802  out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
803  else
804  out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
805  }
806  else if (transpose_matrix == true && nn % 2 == 1)
807  {
808  Number r0;
809  if (m_half > 0)
810  {
811  r0 = matrix[n_half] * xp[0];
812  for (int ind = 1; ind < m_half; ++ind)
813  r0 += matrix[ind * offset + n_half] * xp[ind];
814  }
815  else
816  r0 = Number();
817  if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
818  r0 += matrix[m_half * offset + n_half] * xmid;
819 
820  if (add)
821  out[stride_out * n_half] += r0;
822  else
823  out[stride_out * n_half] = r0;
824  }
825  else if (transpose_matrix == false && nn % 2 == 1)
826  {
827  Number r0;
828  if (m_half > 0)
829  {
830  if (quantity == EvaluatorQuantity::gradient)
831  {
832  r0 = matrix[n_half * offset] * xm[0];
833  for (int ind = 1; ind < m_half; ++ind)
834  r0 += matrix[n_half * offset + ind] * xm[ind];
835  }
836  else
837  {
838  r0 = matrix[n_half * offset] * xp[0];
839  for (int ind = 1; ind < m_half; ++ind)
840  r0 += matrix[n_half * offset + ind] * xp[ind];
841  }
842  }
843  else
844  r0 = Number();
845 
846  if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
847  r0 += matrix[n_half * offset + m_half] * xmid;
848 
849  if (add)
850  out[stride_out * n_half] += r0;
851  else
852  out[stride_out * n_half] = r0;
853  }
854  }
855 
856 
857 
862  template <EvaluatorVariant variant,
863  EvaluatorQuantity quantity,
864  bool transpose_matrix,
865  bool add,
866  typename Number,
867  typename Number2>
868  std::enable_if_t<(variant == evaluate_evenodd), void>
870  const Number *in,
871  Number *out,
872  int n_rows,
873  int n_columns,
874  int stride_in,
875  int stride_out)
876  {
878  quantity,
879  0,
880  0,
881  0,
882  0,
883  transpose_matrix,
884  add>(
885  matrix, in, out, n_rows, n_columns, stride_in, stride_out);
886  }
887 
888 
889 
905  template <EvaluatorVariant variant,
906  EvaluatorQuantity quantity,
907  int n_rows,
908  int n_columns,
909  int stride_in,
910  int stride_out,
911  bool transpose_matrix,
912  bool add,
913  typename Number,
914  typename Number2>
915  std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void>
917  const Number *in,
918  Number *out)
919  {
920  static_assert(n_rows > 0 && n_columns > 0,
921  "Specialization requires n_rows, n_columns > 0");
922 
923  constexpr bool evaluate_antisymmetric =
924  (quantity == EvaluatorQuantity::gradient);
925 
926  constexpr int mm = transpose_matrix ? n_rows : n_columns,
927  nn = transpose_matrix ? n_columns : n_rows;
928  constexpr int n_half = nn / 2;
929  constexpr int m_half = mm / 2;
930 
931  if (transpose_matrix)
932  {
933  std::array<Number, mm> x;
934  for (unsigned int i = 0; i < mm; ++i)
935  x[i] = in[stride_in * i];
936  for (unsigned int col = 0; col < n_half; ++col)
937  {
938  Number r0, r1;
939  if (m_half > 0)
940  {
941  r0 = matrix[col] * x[0];
942  r1 = matrix[col + n_columns] * x[1];
943  for (unsigned int ind = 1; ind < m_half; ++ind)
944  {
945  r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
946  r1 +=
947  matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
948  }
949  }
950  else
951  r0 = r1 = Number();
952  if (mm % 2 == 1)
953  r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
954  if (add)
955  {
956  out[stride_out * col] += r0 + r1;
957  if (evaluate_antisymmetric)
958  out[stride_out * (nn - 1 - col)] += r1 - r0;
959  else
960  out[stride_out * (nn - 1 - col)] += r0 - r1;
961  }
962  else
963  {
964  out[stride_out * col] = r0 + r1;
965  if (evaluate_antisymmetric)
966  out[stride_out * (nn - 1 - col)] = r1 - r0;
967  else
968  out[stride_out * (nn - 1 - col)] = r0 - r1;
969  }
970  }
971  if (nn % 2 == 1)
972  {
973  Number r0;
974  const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
975  if (m_half > 0)
976  {
977  r0 = matrix[n_half + shift * n_columns] * x[shift];
978  for (unsigned int ind = 1; ind < m_half; ++ind)
979  r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
980  x[2 * ind + shift];
981  }
982  else
983  r0 = 0;
984  if (!evaluate_antisymmetric && mm % 2 == 1)
985  r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
986  if (add)
987  out[stride_out * n_half] += r0;
988  else
989  out[stride_out * n_half] = r0;
990  }
991  }
992  else
993  {
994  std::array<Number, m_half + 1> xp, xm;
995  for (int i = 0; i < m_half; ++i)
996  if (!evaluate_antisymmetric)
997  {
998  xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
999  xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1000  }
1001  else
1002  {
1003  xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1004  xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1005  }
1006  if (mm % 2 == 1)
1007  xp[m_half] = in[stride_in * m_half];
1008  for (unsigned int col = 0; col < n_half; ++col)
1009  {
1010  Number r0, r1;
1011  if (m_half > 0)
1012  {
1013  r0 = matrix[2 * col * n_columns] * xp[0];
1014  r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1015  for (unsigned int ind = 1; ind < m_half; ++ind)
1016  {
1017  r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1018  r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1019  }
1020  }
1021  else
1022  r0 = r1 = Number();
1023  if (mm % 2 == 1)
1024  {
1025  if (evaluate_antisymmetric)
1026  r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1027  else
1028  r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1029  }
1030  if (add)
1031  {
1032  out[stride_out * (2 * col)] += r0;
1033  out[stride_out * (2 * col + 1)] += r1;
1034  }
1035  else
1036  {
1037  out[stride_out * (2 * col)] = r0;
1038  out[stride_out * (2 * col + 1)] = r1;
1039  }
1040  }
1041  if (nn % 2 == 1)
1042  {
1043  Number r0;
1044  if (m_half > 0)
1045  {
1046  r0 = matrix[(nn - 1) * n_columns] * xp[0];
1047  for (unsigned int ind = 1; ind < m_half; ++ind)
1048  r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1049  }
1050  else
1051  r0 = Number();
1052  if (mm % 2 == 1 && !evaluate_antisymmetric)
1053  r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1054  if (add)
1055  out[stride_out * (nn - 1)] += r0;
1056  else
1057  out[stride_out * (nn - 1)] = r0;
1058  }
1059  }
1060  }
1061 
1062 
1063 
1086  template <EvaluatorVariant variant,
1087  int dim,
1088  int n_rows,
1089  int n_columns,
1090  typename Number,
1091  typename Number2 = Number>
1093  {
1094  static constexpr unsigned int n_rows_of_product =
1095  Utilities::pow(n_rows, dim);
1096  static constexpr unsigned int n_columns_of_product =
1097  Utilities::pow(n_columns, dim);
1098 
1104  : shape_values(nullptr)
1105  , shape_gradients(nullptr)
1106  , shape_hessians(nullptr)
1107  {}
1108 
1115  const unsigned int = 0,
1116  const unsigned int = 0)
1120  {
1121  if (variant == evaluate_evenodd)
1122  {
1123  if (!shape_values.empty())
1125  n_rows * ((n_columns + 1) / 2));
1126  if (!shape_gradients.empty())
1128  n_rows * ((n_columns + 1) / 2));
1129  if (!shape_hessians.empty())
1131  n_rows * ((n_columns + 1) / 2));
1132  }
1133  else
1134  {
1135  Assert(shape_values.empty() ||
1136  shape_values.size() == n_rows * n_columns,
1137  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1138  Assert(shape_gradients.empty() ||
1139  shape_gradients.size() == n_rows * n_columns,
1141  n_rows * n_columns));
1142  Assert(shape_hessians.empty() ||
1143  shape_hessians.size() == n_rows * n_columns,
1145  n_rows * n_columns));
1146  }
1147  }
1148 
1153  const Number2 *shape_gradients,
1154  const Number2 *shape_hessians,
1155  const unsigned int dummy1 = 0,
1156  const unsigned int dummy2 = 0)
1160  {
1161  (void)dummy1;
1162  (void)dummy2;
1163  }
1164 
1190  template <int direction, bool contract_over_rows, bool add, int stride = 1>
1191  void
1192  values(const Number in[], Number out[]) const
1193  {
1194  constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1195  apply<direction, contract_over_rows, add, false, value_type, stride>(
1196  shape_values, in, out);
1197  }
1198 
1204  template <int direction, bool contract_over_rows, bool add, int stride = 1>
1205  void
1206  gradients(const Number in[], Number out[]) const
1207  {
1208  constexpr EvaluatorQuantity gradient_type =
1211  apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1212  shape_gradients, in, out);
1213  }
1214 
1220  template <int direction, bool contract_over_rows, bool add>
1221  void
1222  hessians(const Number in[], Number out[]) const
1223  {
1224  constexpr EvaluatorQuantity hessian_type =
1225  (((variant == evaluate_general) |
1226  (variant == evaluate_symmetric_hierarchical)) ?
1229  apply<direction, contract_over_rows, add, false, hessian_type>(
1230  shape_hessians, in, out);
1231  }
1232 
1240  template <int direction, bool contract_over_rows, bool add>
1241  void
1242  values_one_line(const Number in[], Number out[]) const
1243  {
1244  Assert(shape_values != nullptr, ExcNotInitialized());
1245  apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1246  shape_values, in, out);
1247  }
1248 
1256  template <int direction, bool contract_over_rows, bool add>
1257  void
1258  gradients_one_line(const Number in[], Number out[]) const
1259  {
1260  Assert(shape_gradients != nullptr, ExcNotInitialized());
1261  constexpr EvaluatorQuantity gradient_type =
1264  apply<direction, contract_over_rows, add, true, gradient_type>(
1265  shape_gradients, in, out);
1266  }
1267 
1275  template <int direction, bool contract_over_rows, bool add>
1276  void
1277  hessians_one_line(const Number in[], Number out[]) const
1278  {
1279  Assert(shape_hessians != nullptr, ExcNotInitialized());
1280  constexpr EvaluatorQuantity hessian_type =
1281  (((variant == evaluate_general) |
1282  (variant == evaluate_symmetric_hierarchical)) ?
1285  apply<direction, contract_over_rows, add, true, hessian_type>(
1286  shape_hessians, in, out);
1287  }
1288 
1325  template <int direction,
1326  bool contract_over_rows,
1327  bool add,
1328  bool one_line = false,
1330  int stride = 1>
1331  static void
1332  apply(const Number2 *DEAL_II_RESTRICT shape_data,
1333  const Number *in,
1334  Number *out);
1335 
1336  private:
1337  const Number2 *shape_values;
1338  const Number2 *shape_gradients;
1339  const Number2 *shape_hessians;
1340  };
1341 
1342 
1343 
1344  template <EvaluatorVariant variant,
1345  int dim,
1346  int n_rows,
1347  int n_columns,
1348  typename Number,
1349  typename Number2>
1350  template <int direction,
1351  bool contract_over_rows,
1352  bool add,
1353  bool one_line,
1354  EvaluatorQuantity quantity,
1355  int stride>
1356  inline void
1358  apply(const Number2 *DEAL_II_RESTRICT shape_data,
1359  const Number *in,
1360  Number *out)
1361  {
1362  static_assert(one_line == false || direction == dim - 1,
1363  "Single-line evaluation only works for direction=dim-1.");
1364  Assert(shape_data != nullptr,
1365  ExcMessage(
1366  "The given array shape_data must not be the null pointer!"));
1367  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1368  in != out,
1369  ExcMessage("In-place operation only supported for "
1370  "n_rows==n_columns or single-line interpolation"));
1371  AssertIndexRange(direction, dim);
1372  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1373  nn = contract_over_rows ? n_columns : n_rows;
1374 
1375  constexpr int stride_operation = Utilities::pow(n_columns, direction);
1376  constexpr int n_blocks1 = one_line ? 1 : stride_operation;
1377  constexpr int n_blocks2 =
1378  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1379 
1380  constexpr int stride_in = !contract_over_rows ? stride : 1;
1381  constexpr int stride_out = contract_over_rows ? stride : 1;
1382  for (int i2 = 0; i2 < n_blocks2; ++i2)
1383  {
1384  for (int i1 = 0; i1 < n_blocks1; ++i1)
1385  {
1387  quantity,
1388  n_rows,
1389  n_columns,
1390  stride_operation * stride_in,
1391  stride_operation * stride_out,
1392  contract_over_rows,
1393  add>(shape_data, in, out);
1394 
1395  if (one_line == false)
1396  {
1397  in += stride_in;
1398  out += stride_out;
1399  }
1400  }
1401  if (one_line == false)
1402  {
1403  in += stride_operation * (mm - 1) * stride_in;
1404  out += stride_operation * (nn - 1) * stride_out;
1405  }
1406  }
1407  }
1408 
1409 
1410 
1424  template <EvaluatorVariant variant,
1425  int dim,
1426  typename Number,
1427  typename Number2>
1428  struct EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>
1429  {
1430  static constexpr unsigned int n_rows_of_product =
1432  static constexpr unsigned int n_columns_of_product =
1434 
1440  : shape_values(nullptr)
1441  , shape_gradients(nullptr)
1442  , shape_hessians(nullptr)
1443  , n_rows(numbers::invalid_unsigned_int)
1444  , n_columns(numbers::invalid_unsigned_int)
1445  {}
1446 
1453  const unsigned int n_rows = 0,
1454  const unsigned int n_columns = 0)
1458  , n_rows(n_rows)
1459  , n_columns(n_columns)
1460  {
1461  if (variant == evaluate_evenodd)
1462  {
1463  if (!shape_values.empty())
1465  n_rows * ((n_columns + 1) / 2));
1466  if (!shape_gradients.empty())
1468  n_rows * ((n_columns + 1) / 2));
1469  if (!shape_hessians.empty())
1471  n_rows * ((n_columns + 1) / 2));
1472  }
1473  else
1474  {
1475  Assert(shape_values.empty() ||
1476  shape_values.size() == n_rows * n_columns,
1477  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1478  Assert(shape_gradients.empty() ||
1479  shape_gradients.size() == n_rows * n_columns,
1481  n_rows * n_columns));
1482  Assert(shape_hessians.empty() ||
1483  shape_hessians.size() == n_rows * n_columns,
1485  n_rows * n_columns));
1486  }
1487  }
1488 
1493  const Number2 *shape_gradients,
1494  const Number2 *shape_hessians,
1495  const unsigned int n_rows = 0,
1496  const unsigned int n_columns = 0)
1500  , n_rows(n_rows)
1501  , n_columns(n_columns)
1502  {}
1503 
1504  template <int direction, bool contract_over_rows, bool add, int stride = 1>
1505  void
1506  values(const Number *in, Number *out) const
1507  {
1508  constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1509  apply<direction, contract_over_rows, add, false, value_type, stride>(
1510  shape_values, in, out);
1511  }
1512 
1513  template <int direction, bool contract_over_rows, bool add, int stride = 1>
1514  void
1515  gradients(const Number *in, Number *out) const
1516  {
1517  constexpr EvaluatorQuantity gradient_type =
1520  apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1521  shape_gradients, in, out);
1522  }
1523 
1524  template <int direction, bool contract_over_rows, bool add>
1525  void
1526  hessians(const Number *in, Number *out) const
1527  {
1528  constexpr EvaluatorQuantity hessian_type =
1531  apply<direction, contract_over_rows, add, false, hessian_type>(
1532  shape_hessians, in, out);
1533  }
1534 
1535  template <int direction, bool contract_over_rows, bool add>
1536  void
1537  values_one_line(const Number in[], Number out[]) const
1538  {
1539  Assert(shape_values != nullptr, ExcNotInitialized());
1540  apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1541  shape_values, in, out);
1542  }
1543 
1544  template <int direction, bool contract_over_rows, bool add>
1545  void
1546  gradients_one_line(const Number in[], Number out[]) const
1547  {
1548  Assert(shape_gradients != nullptr, ExcNotInitialized());
1549  constexpr EvaluatorQuantity gradient_type =
1552  apply<direction, contract_over_rows, add, true, gradient_type>(
1553  shape_gradients, in, out);
1554  }
1555 
1556  template <int direction, bool contract_over_rows, bool add>
1557  void
1558  hessians_one_line(const Number in[], Number out[]) const
1559  {
1560  Assert(shape_hessians != nullptr, ExcNotInitialized());
1561  constexpr EvaluatorQuantity hessian_type =
1564  apply<direction, contract_over_rows, add, true, hessian_type>(
1565  shape_hessians, in, out);
1566  }
1567 
1568  template <int direction,
1569  bool contract_over_rows,
1570  bool add,
1571  bool one_line = false,
1573  int stride = 1>
1574  void
1575  apply(const Number2 *DEAL_II_RESTRICT shape_data,
1576  const Number *in,
1577  Number *out) const;
1578 
1579  const Number2 *shape_values;
1580  const Number2 *shape_gradients;
1581  const Number2 *shape_hessians;
1582  const unsigned int n_rows;
1583  const unsigned int n_columns;
1584  };
1585 
1586 
1587 
1588  template <EvaluatorVariant variant,
1589  int dim,
1590  typename Number,
1591  typename Number2>
1592  template <int direction,
1593  bool contract_over_rows,
1594  bool add,
1595  bool one_line,
1596  EvaluatorQuantity quantity,
1597  int stride>
1598  inline void
1600  const Number2 *DEAL_II_RESTRICT shape_data,
1601  const Number *in,
1602  Number *out) const
1603  {
1604  static_assert(one_line == false || direction == dim - 1,
1605  "Single-line evaluation only works for direction=dim-1.");
1606  Assert(shape_data != nullptr,
1607  ExcMessage(
1608  "The given array shape_data must not be the null pointer!"));
1609  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1610  in != out,
1611  ExcMessage("In-place operation only supported for "
1612  "n_rows==n_columns or single-line interpolation"));
1613  AssertIndexRange(direction, dim);
1614  const int mm = contract_over_rows ? n_rows : n_columns,
1615  nn = contract_over_rows ? n_columns : n_rows;
1616 
1617  const int stride_operation =
1618  direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
1619  const int n_blocks1 = one_line ? 1 : stride_operation;
1620  const int n_blocks2 = direction >= dim - 1 ?
1621  1 :
1622  Utilities::fixed_power<dim - direction - 1>(n_rows);
1623  Assert(n_rows <= 128, ExcNotImplemented());
1624 
1625  constexpr int stride_in = !contract_over_rows ? stride : 1;
1626  constexpr int stride_out = contract_over_rows ? stride : 1;
1627  for (int i2 = 0; i2 < n_blocks2; ++i2)
1628  {
1629  for (int i1 = 0; i1 < n_blocks1; ++i1)
1630  {
1631  // the empty template case can only run the general evaluator or
1632  // evenodd
1633  constexpr EvaluatorVariant restricted_variant =
1635  apply_matrix_vector_product<restricted_variant,
1636  quantity,
1637  contract_over_rows,
1638  add>(shape_data,
1639  in,
1640  out,
1641  n_rows,
1642  n_columns,
1643  stride_operation * stride_in,
1644  stride_operation * stride_out);
1645 
1646  if (one_line == false)
1647  {
1648  in += stride_in;
1649  out += stride_out;
1650  }
1651  }
1652  if (one_line == false)
1653  {
1654  in += stride_operation * (mm - 1) * stride_in;
1655  out += stride_operation * (nn - 1) * stride_out;
1656  }
1657  }
1658  }
1659 
1660 
1661 
1662  template <int dim,
1663  int fe_degree,
1664  int n_q_points_1d,
1665  bool contract_over_rows,
1666  bool symmetric_evaluate = true>
1668  {
1669  template <int direction,
1670  int stride = 1,
1671  typename Number = double,
1672  typename Number2 = double>
1673  static void
1675  const Number *in,
1676  Number *out,
1677  const bool add_into_result = false,
1678  const int subface_index_1d = 0)
1679  {
1680  AssertIndexRange(direction, dim);
1681  AssertDimension(fe_degree, data.fe_degree);
1682  AssertDimension(n_q_points_1d, data.n_q_points_1d);
1683  constexpr int n_rows = fe_degree + 1;
1684  constexpr int n_columns = n_q_points_1d;
1685  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1686  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1687  const Number2 *shape_data =
1688  symmetric_evaluate ?
1689  data.shape_values_eo.data() :
1690  data.values_within_subface[subface_index_1d].data();
1691  Assert(shape_data != nullptr, ExcNotInitialized());
1692  Assert(contract_over_rows == false || !add_into_result,
1693  ExcMessage("Cannot add into result if contract_over_rows = true"));
1694 
1695  constexpr int n_blocks1 = Utilities::pow(fe_degree, direction);
1696  constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1);
1697  constexpr int stride_in = contract_over_rows ? 1 : stride;
1698  constexpr int stride_out = contract_over_rows ? stride : 1;
1699  constexpr EvaluatorVariant variant =
1700  symmetric_evaluate ? evaluate_evenodd : evaluate_general;
1701 
1702  for (int i2 = 0; i2 < n_blocks2; ++i2)
1703  {
1704  for (int i1 = 0; i1 < n_blocks1; ++i1)
1705  {
1706  if (contract_over_rows == false && add_into_result)
1709  n_rows,
1710  n_columns,
1711  n_blocks1 * stride_in,
1712  n_blocks1 * stride_out,
1713  contract_over_rows,
1714  true>(shape_data, in, out);
1715  else
1718  n_rows,
1719  n_columns,
1720  n_blocks1 * stride_in,
1721  n_blocks1 * stride_out,
1722  contract_over_rows,
1723  false>(shape_data, in, out);
1724 
1725  in += stride_in;
1726  out += stride_out;
1727  }
1728  in += n_blocks1 * (mm - 1) * stride_in;
1729  out += n_blocks1 * (nn - 1) * stride_out;
1730  }
1731  }
1732 
1733  template <int direction,
1734  int normal_direction,
1735  int stride = 1,
1736  typename Number = double,
1737  typename Number2 = double>
1738  static void
1740  const Number *in,
1741  Number *out,
1742  const int subface_index_1d = 0)
1743  {
1744  AssertIndexRange(direction, dim);
1745  AssertDimension(fe_degree - 1, data.fe_degree);
1746  AssertDimension(n_q_points_1d, data.n_q_points_1d);
1747  static_assert(direction != normal_direction,
1748  "Cannot interpolate tangentially in normal direction");
1749 
1750  constexpr int n_rows = std::max(fe_degree, 0);
1751  constexpr int n_columns = n_q_points_1d;
1752  const Number2 *shape_data =
1753  symmetric_evaluate ?
1754  data.shape_values_eo.data() :
1755  data.values_within_subface[subface_index_1d].data();
1756  Assert(shape_data != nullptr, ExcNotInitialized());
1757 
1758  constexpr int n_blocks1 =
1759  (direction > normal_direction) ?
1760  Utilities::pow(n_q_points_1d, direction) :
1761  (direction > 0 ?
1762  (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
1763  1);
1764  constexpr int n_blocks2 =
1765  (direction > normal_direction) ?
1766  Utilities::pow(fe_degree, dim - 1 - direction) :
1767  ((direction + 1 < dim) ?
1768  (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
1769  1);
1770 
1771  constexpr EvaluatorVariant variant =
1772  symmetric_evaluate ? evaluate_evenodd : evaluate_general;
1773 
1774  // Since we may perform an in-place interpolation, we must run the step
1775  // expanding the size of the basis backward ('contract_over_rows' aka
1776  // 'evaluate' case), so shift the pointers and decrement during the loop
1777  if (contract_over_rows)
1778  {
1779  in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
1780  out +=
1781  stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
1782  for (int i2 = 0; i2 < n_blocks2; ++i2)
1783  {
1784  for (int i1 = 0; i1 < n_blocks1; ++i1)
1785  {
1788  n_rows,
1789  n_columns,
1790  n_blocks1,
1791  n_blocks1 * stride,
1792  true,
1793  false>(shape_data, in, out);
1794 
1795  --in;
1796  out -= stride;
1797  }
1798  in -= n_blocks1 * (n_rows - 1);
1799  out -= n_blocks1 * (n_columns - 1) * stride;
1800  }
1801  }
1802  else
1803  {
1804  for (int i2 = 0; i2 < n_blocks2; ++i2)
1805  {
1806  for (int i1 = 0; i1 < n_blocks1; ++i1)
1807  {
1810  n_rows,
1811  n_columns,
1812  n_blocks1 * stride,
1813  n_blocks1,
1814  false,
1815  false>(shape_data, in, out);
1816 
1817  in += stride;
1818  ++out;
1819  }
1820  in += n_blocks1 * (n_columns - 1) * stride;
1821  out += n_blocks1 * (n_rows - 1);
1822  }
1823  }
1824  }
1825  };
1826 
1827 
1828 
1875  template <int n_rows_template,
1876  int stride_template,
1877  bool contract_onto_face,
1878  bool add,
1879  int max_derivative,
1880  typename Number,
1881  typename Number2>
1882  inline std::enable_if_t<contract_onto_face, void>
1883  interpolate_to_face(const Number2 *shape_values,
1884  const std::array<int, 2> &n_blocks,
1885  const std::array<int, 2> &steps,
1886  const Number *input,
1887  Number *DEAL_II_RESTRICT output,
1888  const int n_rows_runtime = 0,
1889  const int stride_runtime = 1)
1890  {
1891  const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
1892  const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
1893 
1894  Number *output1 = output + n_blocks[0] * n_blocks[1];
1895  Number *output2 = output1 + n_blocks[0] * n_blocks[1];
1896  for (int i2 = 0; i2 < n_blocks[1]; ++i2)
1897  {
1898  for (int i1 = 0; i1 < n_blocks[0]; ++i1)
1899  {
1900  Number res0 = shape_values[0] * input[0];
1901  Number res1, res2;
1902  if (max_derivative > 0)
1903  res1 = shape_values[n_rows] * input[0];
1904  if (max_derivative > 1)
1905  res2 = shape_values[2 * n_rows] * input[0];
1906  for (int ind = 1; ind < n_rows; ++ind)
1907  {
1908  res0 += shape_values[ind] * input[stride * ind];
1909  if (max_derivative > 0)
1910  res1 += shape_values[ind + n_rows] * input[stride * ind];
1911  if (max_derivative > 1)
1912  res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
1913  }
1914  if (add)
1915  {
1916  output[i1] += res0;
1917  if (max_derivative > 0)
1918  output1[i1] += res1;
1919  if (max_derivative > 1)
1920  output2[i2] += res2;
1921  }
1922  else
1923  {
1924  output[i1] = res0;
1925  if (max_derivative > 0)
1926  output1[i1] = res1;
1927  if (max_derivative > 1)
1928  output2[i1] = res2;
1929  }
1930  input += steps[0];
1931  }
1932  output += n_blocks[0];
1933  if (max_derivative > 0)
1934  output1 += n_blocks[0];
1935  if (max_derivative > 1)
1936  output2 += n_blocks[0];
1937  input += steps[1];
1938  }
1939  }
1940 
1941 
1942 
1950  constexpr bool
1951  use_collocation_evaluation(const unsigned int fe_degree,
1952  const unsigned int n_q_points_1d)
1953  {
1954  return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
1955  (n_q_points_1d <= 3 * fe_degree / 2 + 1);
1956  }
1957 
1958 
1959 
1965  template <int n_rows_template,
1966  int stride_template,
1967  bool contract_onto_face,
1968  bool add,
1969  int max_derivative,
1970  typename Number,
1971  typename Number2>
1972  inline std::enable_if_t<!contract_onto_face, void>
1973  interpolate_to_face(const Number2 *shape_values,
1974  const std::array<int, 2> &n_blocks,
1975  const std::array<int, 2> &steps,
1976  const Number *input,
1977  Number *DEAL_II_RESTRICT output,
1978  const int n_rows_runtime = 0,
1979  const int stride_runtime = 1)
1980  {
1981  const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
1982  const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
1983 
1984  const Number *input1 = input + n_blocks[0] * n_blocks[1];
1985  const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
1986  for (int i2 = 0; i2 < n_blocks[1]; ++i2)
1987  {
1988  for (int i1 = 0; i1 < n_blocks[0]; ++i1)
1989  {
1990  const Number in = input[i1];
1991  Number in1, in2;
1992  if (max_derivative > 0)
1993  in1 = input1[i1];
1994  if (max_derivative > 1)
1995  in2 = input2[i1];
1996  for (int col = 0; col < n_rows; ++col)
1997  {
1998  Number result =
1999  add ? (output[col * stride] + shape_values[col] * in) :
2000  (shape_values[col] * in);
2001  if (max_derivative > 0)
2002  result += shape_values[col + n_rows] * in1;
2003  if (max_derivative > 1)
2004  result += shape_values[col + 2 * n_rows] * in2;
2005 
2006  output[col * stride] = result;
2007  }
2008  output += steps[0];
2009  }
2010  input += n_blocks[0];
2011  if (max_derivative > 0)
2012  input1 += n_blocks[0];
2013  if (max_derivative > 1)
2014  input2 += n_blocks[0];
2015  output += steps[1];
2016  }
2017  }
2018 
2019 
2020 
2027  template <typename Number, typename Number2>
2029  {
2031  };
2032 
2033  template <int dim, typename Number, typename Number2>
2034  struct ProductTypeNoPoint<Point<dim, Number>, Number2>
2035  {
2037  };
2038 
2039 
2040 
2045  template <int dim, typename Number>
2046  inline void
2048  ::ndarray<Number, 2, dim> *shapes,
2049  const std::vector<Polynomials::Polynomial<double>> &poly,
2050  const Point<dim, Number> &p,
2051  const unsigned int derivative = 1)
2052  {
2053  const int n_shapes = poly.size();
2054 
2055  // Evaluate 1d polynomials and their derivatives
2056  std::array<Number, dim> point;
2057  for (unsigned int d = 0; d < dim; ++d)
2058  point[d] = p[d];
2059  for (int i = 0; i < n_shapes; ++i)
2060  poly[i].values_of_array(point, derivative, shapes[i].data());
2061  }
2062 
2063 
2064 
2068  template <typename Number>
2069  inline void
2071  const std::vector<Polynomials::Polynomial<double>> &,
2072  const Point<0, Number> &,
2073  const unsigned int)
2074  {
2075  Assert(false, ExcInternalError());
2076  }
2077 
2078 
2079 
2083  template <int dim,
2084  int length,
2085  typename Number2,
2086  typename Number,
2087  int n_values = 1,
2088  bool do_renumber = true>
2089  inline
2090 #ifndef DEBUG
2092 #endif
2093  std::array<typename ProductTypeNoPoint<Number, Number2>::type,
2094  2 + n_values>
2095  do_interpolate_xy(const Number *values,
2096  const std::vector<unsigned int> &renumber,
2097  const ::ndarray<Number2, 2, dim> *shapes,
2098  const int n_shapes_runtime,
2099  int &i)
2100  {
2101  static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
2102  static_assert(1 <= n_values && n_values <= 2,
2103  "Only n_values=1,2 implemented");
2104 
2105  const int n_shapes = length > 0 ? length : n_shapes_runtime;
2106 
2107  // If n_values > 1, we want to interpolate from a second array,
2108  // placed in the same array immediately after the main data. This
2109  // is used to interpolate normal derivatives onto faces.
2110  const Number *values_2 =
2111  n_values > 1 ?
2112  values + (length > 0 ? Utilities::pow(length, dim) :
2113  Utilities::fixed_power<dim>(n_shapes_runtime)) :
2114  nullptr;
2115  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2116  std::array<Number3, 2 + n_values> result = {};
2117  for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2118  {
2119  // Interpolation + derivative x direction
2120  std::array<Number3, 1 + n_values> inner_result = {};
2121 
2122  // Distinguish the inner loop based on whether we have a
2123  // renumbering or not
2124  if (do_renumber && !renumber.empty())
2125  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2126  {
2127  // gradient
2128  inner_result[0] += shapes[i0][1][0] * values[renumber[i]];
2129  // values
2130  inner_result[1] += shapes[i0][0][0] * values[renumber[i]];
2131  if (n_values > 1)
2132  inner_result[2] += shapes[i0][0][0] * values_2[renumber[i]];
2133  }
2134  else
2135  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2136  {
2137  // gradient
2138  inner_result[0] += shapes[i0][1][0] * values[i];
2139  // values
2140  inner_result[1] += shapes[i0][0][0] * values[i];
2141  if (n_values > 1)
2142  inner_result[2] += shapes[i0][0][0] * values_2[i];
2143  }
2144 
2145  if (dim > 1)
2146  {
2147  // Interpolation + derivative in y direction
2148  // gradient
2149  result[0] += inner_result[0] * shapes[i1][0][1];
2150  result[1] += inner_result[1] * shapes[i1][1][1];
2151  // values
2152  result[2] += inner_result[1] * shapes[i1][0][1];
2153  if (n_values > 1)
2154  result[3] += inner_result[2] * shapes[i1][0][1];
2155  }
2156  else
2157  {
2158  // gradient
2159  result[0] = inner_result[0];
2160  // values
2161  result[1] = inner_result[1];
2162  if (n_values > 1)
2163  result[2] = inner_result[2];
2164  }
2165  }
2166  return result;
2167  }
2168 
2169 
2170 
2175  template <int dim,
2176  typename Number,
2177  typename Number2,
2178  int n_values = 1,
2179  bool do_renumber = true>
2180  inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
2181  dim + n_values>
2183  const ::ndarray<Number2, 2, dim> *shapes,
2184  const int n_shapes,
2185  const Number *values,
2186  const std::vector<unsigned int> &renumber = {})
2187  {
2188  static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
2189  static_assert(1 <= n_values && n_values <= 2,
2190  "Only n_values=1,2 implemented");
2191 
2192  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2193 
2194  std::array<Number3, dim + n_values> result = {};
2195  if (dim == 0)
2196  {
2197  // We only need the interpolation of the value and normal derivatives on
2198  // faces of a 1d element. As the interpolation is the value at the
2199  // point, simply set the result vector accordingly.
2200  result[0] = values[0];
2201  if (n_values > 1)
2202  result[1] = values[1];
2203  return result;
2204  }
2205 
2206  // Go through the tensor product of shape functions and interpolate
2207  // with optimal algorithm
2208  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2209  {
2210  std::array<Number3, 2 + n_values> inner_result;
2211  // Generate separate code with known loop bounds for the most common
2212  // cases
2213  if (n_shapes == 2)
2214  inner_result =
2215  do_interpolate_xy<dim, 2, Number2, Number, n_values, do_renumber>(
2216  values, renumber, shapes, n_shapes, i);
2217  else if (n_shapes == 3)
2218  inner_result =
2219  do_interpolate_xy<dim, 3, Number2, Number, n_values, do_renumber>(
2220  values, renumber, shapes, n_shapes, i);
2221  else if (n_shapes == 4)
2222  inner_result =
2223  do_interpolate_xy<dim, 4, Number2, Number, n_values, do_renumber>(
2224  values, renumber, shapes, n_shapes, i);
2225  else if (n_shapes == 5)
2226  inner_result =
2227  do_interpolate_xy<dim, 5, Number2, Number, n_values, do_renumber>(
2228  values, renumber, shapes, n_shapes, i);
2229  else if (n_shapes == 6)
2230  inner_result =
2231  do_interpolate_xy<dim, 6, Number2, Number, n_values, do_renumber>(
2232  values, renumber, shapes, n_shapes, i);
2233  else
2234  inner_result =
2235  do_interpolate_xy<dim, -1, Number2, Number, n_values, do_renumber>(
2236  values, renumber, shapes, n_shapes, i);
2237  if (dim == 3)
2238  {
2239  // derivative + interpolation in z direction
2240  // gradient
2241  result[0] += inner_result[0] * shapes[i2][0][2];
2242  result[1] += inner_result[1] * shapes[i2][0][2];
2243  result[2] += inner_result[2] * shapes[i2][1][2];
2244  // values
2245  result[3] += inner_result[2] * shapes[i2][0][2];
2246  if (n_values > 1)
2247  result[4] += inner_result[3] * shapes[i2][0][2];
2248  }
2249  else if (dim == 2)
2250  {
2251  // gradient
2252  result[0] = inner_result[0];
2253  result[1] = inner_result[1];
2254  // values
2255  result[2] = inner_result[2];
2256  if (n_values > 1)
2257  result[3] = inner_result[3];
2258  }
2259  else
2260  {
2261  // gradient
2262  result[0] = inner_result[0];
2263  // values
2264  result[1] = inner_result[1];
2265  if (n_values > 1)
2266  result[2] = inner_result[2];
2267  }
2268  }
2269 
2270  return result;
2271  }
2272 
2273 
2274 
2279  template <int dim,
2280  typename Number,
2281  typename Number2,
2282  int n_values = 1,
2283  int stride = 1>
2284  inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
2285  dim + n_values>
2287  const Number *values,
2288  const Point<dim, Number2> &p)
2289  {
2290  static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
2291  static_assert(1 <= n_values && n_values <= 2,
2292  "Only n_values=1,2 implemented");
2293 
2294  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2295 
2296  static_assert(
2297  n_values == 1 || stride == 1,
2298  "Either n_values or stride has to be one for correct data access!");
2299  // If n_values > 1, we want to interpolate from a second array,
2300  // placed in the same array immediately after the main data. This
2301  // is used to interpolate normal derivatives onto faces.
2302 
2303  std::array<Number3, dim + n_values> result;
2304  if (dim == 0)
2305  {
2306  // we only need the value on faces of a 1d element
2307  result[0] = values[0];
2308  if (n_values > 1)
2309  result[1] = values[1];
2310  }
2311  else if (dim == 1)
2312  {
2313  // gradient
2314  result[0] = Number3(values[stride] - values[0]);
2315  // values
2316  result[1] = Number3(values[0]) + p[0] * result[0];
2317  if (n_values > 1)
2318  result[2] = Number3(values[2]) + p[0] * (values[3] - values[2]);
2319  }
2320  else if (dim == 2)
2321  {
2322  const Number3 val10 = Number3(values[stride] - values[0]);
2323  const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
2324  const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
2325  const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
2326 
2327  // gradient
2328  result[0] = val10 + p[1] * (val32 - val10);
2329  result[1] = tmp1 - tmp0;
2330 
2331  // values
2332  result[2] = tmp0 + p[1] * result[1];
2333 
2334  if (n_values > 1)
2335  {
2336  const Number3 tmp0_2 =
2337  Number3(values[4]) + p[0] * (values[5] - values[4]);
2338  const Number3 tmp1_2 =
2339  Number3(values[6]) + p[0] * (values[7] - values[6]);
2340  result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
2341  }
2342  }
2343  else if (dim == 3)
2344  {
2345  const Number3 val10 = Number3(values[stride] - values[0]);
2346  const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
2347  const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
2348  const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
2349  const Number3 tmp10 = tmp1 - tmp0;
2350  const Number3 tmpy0 = tmp0 + p[1] * tmp10;
2351 
2352  const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
2353  const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
2354  const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54;
2355  const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76;
2356  const Number3 tmp32 = tmp3 - tmp2;
2357  const Number3 tmpy1 = tmp2 + p[1] * tmp32;
2358 
2359  // gradient
2360  result[2] = tmpy1 - tmpy0;
2361  result[1] = tmp10 + p[2] * (tmp32 - tmp10);
2362  const Number3 tmpz0 = val10 + p[1] * (val32 - val10);
2363  result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0);
2364 
2365  // value
2366  result[3] = tmpy0 + p[2] * result[2];
2367  Assert(n_values == 1, ExcNotImplemented());
2368  }
2369 
2370  return result;
2371  }
2372 
2373 
2374 
2409  template <int dim, typename Number, typename Number2>
2410  inline std::pair<
2414  const std::vector<Polynomials::Polynomial<double>> &poly,
2415  const std::vector<Number> &values,
2416  const Point<dim, Number2> &p,
2417  const bool d_linear = false,
2418  const std::vector<unsigned int> &renumber = {})
2419  {
2420  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2421 
2422  std::array<Number3, dim + 1> result;
2423  if (d_linear)
2424  {
2425  result =
2427  }
2428  else
2429  {
2430  AssertIndexRange(poly.size(), 200);
2431  std::array<::ndarray<Number2, 2, dim>, 200> shapes;
2432  compute_values_of_array(shapes.data(), poly, p);
2434  Number,
2435  Number2>(
2436  shapes.data(), poly.size(), values.data(), renumber);
2437  }
2438  return std::make_pair(result[dim],
2440  ArrayView<Number3>(result.data(), dim)));
2441  }
2442 
2443 
2444 
2445  template <int dim,
2446  int length,
2447  typename Number2,
2448  typename Number,
2449  bool do_renumber = true>
2450  inline
2451 #ifndef DEBUG
2453 #endif
2456  const std::vector<unsigned int> &renumber,
2457  const ::ndarray<Number2, 2, dim> *shapes,
2458  const int n_shapes_runtime,
2459  int &i)
2460  {
2461  const int n_shapes = length > 0 ? length : n_shapes_runtime;
2462  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2463  Number3 result = {};
2464  for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2465  {
2466  // Interpolation x direction
2467  Number3 value = {};
2468 
2469  // Distinguish the inner loop based on whether we have a
2470  // renumbering or not
2471  if (do_renumber && !renumber.empty())
2472  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2473  value += shapes[i0][0][0] * values[renumber[i]];
2474  else
2475  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2476  value += shapes[i0][0][0] * values[i];
2477 
2478  if (dim > 1)
2479  result += value * shapes[i1][0][1];
2480  else
2481  result = value;
2482  }
2483  return result;
2484  }
2485 
2486 
2487 
2488  template <int dim, typename Number, typename Number2, bool do_renumber = true>
2491  const ::ndarray<Number2, 2, dim> *shapes,
2492  const int n_shapes,
2493  const Number *values,
2494  const std::vector<unsigned int> &renumber = {})
2495  {
2496  static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
2497 
2498  // we only need the value on faces of a 1d element
2499  if (dim == 0)
2500  {
2501  return values[0];
2502  }
2503 
2504  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2505 
2506  // Go through the tensor product of shape functions and interpolate
2507  // with optimal algorithm
2508  Number3 result = {};
2509  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2510  {
2511  Number3 inner_result;
2512  // Generate separate code with known loop bounds for the most common
2513  // cases
2514  if (n_shapes == 2)
2515  inner_result =
2516  do_interpolate_xy_value<dim, 2, Number2, Number, do_renumber>(
2517  values, renumber, shapes, n_shapes, i);
2518  else if (n_shapes == 3)
2519  inner_result =
2520  do_interpolate_xy_value<dim, 3, Number2, Number, do_renumber>(
2521  values, renumber, shapes, n_shapes, i);
2522  else if (n_shapes == 4)
2523  inner_result =
2524  do_interpolate_xy_value<dim, 4, Number2, Number, do_renumber>(
2525  values, renumber, shapes, n_shapes, i);
2526  else if (n_shapes == 5)
2527  inner_result =
2528  do_interpolate_xy_value<dim, 5, Number2, Number, do_renumber>(
2529  values, renumber, shapes, n_shapes, i);
2530  else if (n_shapes == 6)
2531  inner_result =
2532  do_interpolate_xy_value<dim, 6, Number2, Number, do_renumber>(
2533  values, renumber, shapes, n_shapes, i);
2534  else
2535  inner_result =
2536  do_interpolate_xy_value<dim, -1, Number2, Number, do_renumber>(
2537  values, renumber, shapes, n_shapes, i);
2538  if (dim == 3)
2539  {
2540  // Interpolation + derivative in z direction
2541  result += inner_result * shapes[i2][0][2];
2542  }
2543  else
2544  {
2545  result = inner_result;
2546  }
2547  }
2548 
2549  return result;
2550  }
2551 
2552 
2553 
2554  template <int dim, typename Number, typename Number2, int stride = 1>
2557  const Point<dim, Number2> &p)
2558  {
2559  static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
2560 
2561  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2562 
2563  if (dim == 0)
2564  {
2565  // we only need the value on faces of a 1d element
2566  return values[0];
2567  }
2568  else if (dim == 1)
2569  {
2570  return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]);
2571  }
2572  else if (dim == 2)
2573  {
2574  const Number3 val10 = Number3(values[stride] - values[0]);
2575  const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
2576  const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
2577  const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
2578  return tmp0 + p[1] * (tmp1 - tmp0);
2579  }
2580  else if (dim == 3)
2581  {
2582  const Number3 val10 = Number3(values[stride] - values[0]);
2583  const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
2584  const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
2585  const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
2586  const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
2587 
2588  const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
2589  const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
2590  const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54;
2591  const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76;
2592  const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
2593 
2594  return tmpy0 + p[2] * (tmpy1 - tmpy0);
2595  }
2596 
2597  // work around a compile error: missing return statement
2598  return Number3();
2599  }
2600 
2601 
2602 
2603  template <int dim, typename Number, typename Number2>
2606  const std::vector<Polynomials::Polynomial<double>> &poly,
2607  const std::vector<Number> &values,
2608  const Point<dim, Number2> &p,
2609  const bool d_linear = false,
2610  const std::vector<unsigned int> &renumber = {})
2611  {
2613  if (d_linear)
2614  {
2615  result = evaluate_tensor_product_value_linear(values.data(), p);
2616  }
2617  else
2618  {
2619  AssertIndexRange(poly.size(), 200);
2620  std::array<::ndarray<Number2, 2, dim>, 200> shapes;
2621  const int n_shapes = poly.size();
2622  std::array<Number2, dim> point;
2623  for (unsigned int d = 0; d < dim; ++d)
2624  point[d] = p[d];
2625  for (int i = 0; i < n_shapes; ++i)
2626  poly[i].values_of_array(point, 0, shapes[i].data());
2627  result = evaluate_tensor_product_value_shapes<dim, Number, Number2>(
2628  shapes.data(), n_shapes, values.data(), renumber);
2629  }
2630  return result;
2631  }
2632 
2633 
2634 
2639  template <int derivative_order, typename Number, typename Number2>
2642  const std::vector<Polynomials::Polynomial<double>> &poly,
2643  const std::vector<Number> &values,
2644  const Point<1, Number2> &p,
2645  const std::vector<unsigned int> &renumber = {})
2646  {
2647  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2648 
2649  const int n_shapes = poly.size();
2650  AssertDimension(n_shapes, values.size());
2651  Assert(renumber.empty() || renumber.size() == values.size(),
2652  ExcDimensionMismatch(renumber.size(), values.size()));
2653 
2654  std::array<Number2, derivative_order + 1> shapes;
2655  Tensor<1, 1, Number3> result;
2656  if (renumber.empty())
2657  for (int i = 0; i < n_shapes; ++i)
2658  {
2659  poly[i].value(p[0], derivative_order, shapes.data());
2660  result[0] += shapes[derivative_order] * values[i];
2661  }
2662  else
2663  for (int i = 0; i < n_shapes; ++i)
2664  {
2665  poly[i].value(p[0], derivative_order, shapes.data());
2666  result[0] += shapes[derivative_order] * values[renumber[i]];
2667  }
2668  return result;
2669  }
2670 
2671 
2672 
2677  template <int derivative_order, typename Number, typename Number2>
2678  inline Tensor<1,
2679  derivative_order + 1,
2682  const std::vector<Polynomials::Polynomial<double>> &poly,
2683  const std::vector<Number> &values,
2684  const Point<2, Number2> &p,
2685  const std::vector<unsigned int> &renumber = {})
2686  {
2687  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2688  constexpr int dim = 2;
2689 
2690  const int n_shapes = poly.size();
2691  AssertDimension(Utilities::pow(n_shapes, 2), values.size());
2692  Assert(renumber.empty() || renumber.size() == values.size(),
2693  ExcDimensionMismatch(renumber.size(), values.size()));
2694 
2695  AssertIndexRange(n_shapes, 100);
2697  // Evaluate 1d polynomials and their derivatives
2698  std::array<Number2, dim> point;
2699  for (unsigned int d = 0; d < dim; ++d)
2700  point[d] = p[d];
2701  for (int i = 0; i < n_shapes; ++i)
2702  poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
2703 
2705  for (int i1 = 0, i = 0; i1 < n_shapes; ++i1)
2706  {
2708  if (renumber.empty())
2709  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2710  for (unsigned int d = 0; d <= derivative_order; ++d)
2711  result_x[d] += shapes[i0][d][0] * values[i];
2712  else
2713  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2714  for (unsigned int d = 0; d <= derivative_order; ++d)
2715  result_x[d] += shapes[i0][d][0] * values[renumber[i]];
2716 
2717  for (unsigned int d = 0; d <= derivative_order; ++d)
2718  result[d] += shapes[i1][d][1] * result_x[derivative_order - d];
2719  }
2720  return result;
2721  }
2722 
2723 
2724 
2729  template <int derivative_order, typename Number, typename Number2>
2730  inline Tensor<1,
2731  ((derivative_order + 1) * (derivative_order + 2)) / 2,
2734  const std::vector<Polynomials::Polynomial<double>> &poly,
2735  const std::vector<Number> &values,
2736  const Point<3, Number2> &p,
2737  const std::vector<unsigned int> &renumber = {})
2738  {
2739  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2740  constexpr int dim = 3;
2741  constexpr int n_derivatives =
2742  ((derivative_order + 1) * (derivative_order + 2)) / 2;
2743 
2744  const int n_shapes = poly.size();
2745  AssertDimension(Utilities::pow(n_shapes, 3), values.size());
2746  Assert(renumber.empty() || renumber.size() == values.size(),
2747  ExcDimensionMismatch(renumber.size(), values.size()));
2748 
2749  AssertIndexRange(n_shapes, 100);
2751  // Evaluate 1d polynomials and their derivatives
2752  std::array<Number2, dim> point;
2753  for (unsigned int d = 0; d < dim; ++d)
2754  point[d] = p[d];
2755  for (int i = 0; i < n_shapes; ++i)
2756  poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
2757 
2759  for (int i2 = 0, i = 0; i2 < n_shapes; ++i2)
2760  {
2762  for (int i1 = 0; i1 < n_shapes; ++i1)
2763  {
2764  // apply x derivatives
2766  if (renumber.empty())
2767  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2768  for (unsigned int d = 0; d <= derivative_order; ++d)
2769  result_x[d] += shapes[i0][d][0] * values[i];
2770  else
2771  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2772  for (unsigned int d = 0; d <= derivative_order; ++d)
2773  result_x[d] += shapes[i0][d][0] * values[renumber[i]];
2774 
2775  // multiply by y derivatives, sorting them in upper triangular
2776  // matrix, starting with highest global derivative order,
2777  // decreasing the combined order of xy derivatives by one in each
2778  // row, to be combined with z derivatives in the next step
2779  for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
2780  for (unsigned int e = d; e <= derivative_order; ++e, ++c)
2781  result_xy[c] +=
2782  shapes[i1][e - d][1] * result_x[derivative_order - e];
2783  }
2784 
2785  // multiply by z derivatives, starting with highest x derivative
2786  for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
2787  for (unsigned int e = d; e <= derivative_order; ++e, ++c)
2788  result[c] += shapes[i2][d][2] * result_xy[c];
2789  }
2790  return result;
2791  }
2792 
2793 
2794 
2795  template <int dim, typename Number, typename Number2>
2798  const std::vector<Polynomials::Polynomial<double>> &poly,
2799  const std::vector<Number> &values,
2800  const Point<dim, Number2> &p,
2801  const std::vector<unsigned int> &renumber = {})
2802  {
2803  static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
2804 
2805  const auto hessian =
2806  evaluate_tensor_product_higher_derivatives<2>(poly, values, p, renumber);
2807 
2808  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2810  if (dim == 1)
2811  result[0][0] = hessian[0];
2812  else if (dim >= 2)
2813  {
2814  // derivatives in Hessian are xx, xy, yy, xz, yz, zz, so must re-order
2815  // them for 3D
2816  for (unsigned int d = 0, c = 0; d < 2; ++d)
2817  for (unsigned int e = d; e < 2; ++e, ++c)
2818  result[d][e] = hessian[c];
2819  if (dim == 3)
2820  {
2821  for (unsigned int d = 0; d < 2; ++d)
2822  result[d][2] = hessian[3 + d];
2823  result[2][2] = hessian[5];
2824  }
2825  }
2826 
2827  return result;
2828  }
2829 
2830 
2831 
2835  template <int dim,
2836  int length,
2837  typename Number2,
2838  typename Number,
2839  bool add,
2840  int n_values = 1>
2841  inline
2842 #ifndef DEBUG
2844 #endif
2845  void
2847  Number2 *values,
2848  const ::ndarray<Number, 2, dim> *shapes,
2849  const std::array<Number2, 2 + n_values> &test_grads_value,
2850  const int n_shapes_runtime,
2851  int &i)
2852  {
2853  static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
2854  static_assert(1 <= n_values && n_values <= 2,
2855  "Only n_values=1,2 implemented");
2856 
2857  // Note that 'add' is a template argument, so the compiler will remove
2858  // these checks
2859  if (length > 0)
2860  {
2861  constexpr unsigned int array_size = length > 0 ? length : 1;
2862  std::array<Number, array_size> shape_values_x;
2863  std::array<Number, array_size> shape_derivs_x;
2864  for (unsigned int j = 0; j < array_size; ++j)
2865  {
2866  shape_values_x[j] = shapes[j][0][0];
2867  shape_derivs_x[j] = shapes[j][1][0];
2868  }
2869  for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
2870  {
2871  const Number2 test_value_y =
2872  dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
2873  test_grads_value[1] * shapes[i1][1][1]) :
2874  test_grads_value[2];
2875  const Number2 test_grad_xy =
2876  dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
2877  test_grads_value[0];
2878  Number2 test_value_y_2;
2879  if (n_values > 1)
2880  test_value_y_2 = dim > 1 ?
2881  test_grads_value[3] * shapes[i1][0][1] :
2882  test_grads_value[3];
2883 
2884  Number2 *values_ptr = values + i + i1 * length;
2885  Number2 *values_ptr_2 =
2886  n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr;
2887  for (int i0 = 0; i0 < length; ++i0)
2888  {
2889  if (add)
2890  values_ptr[i0] += shape_values_x[i0] * test_value_y;
2891  else
2892  values_ptr[i0] = shape_values_x[i0] * test_value_y;
2893  values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy;
2894  if (n_values > 1)
2895  {
2896  if (add)
2897  values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2;
2898  else
2899  values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2;
2900  }
2901  }
2902  }
2903  i += (dim > 1 ? length * length : length);
2904  }
2905  else
2906  {
2907  for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
2908  {
2909  const Number2 test_value_y =
2910  dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
2911  test_grads_value[1] * shapes[i1][1][1]) :
2912  test_grads_value[2];
2913  const Number2 test_grad_xy =
2914  dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
2915  test_grads_value[0];
2916  Number2 test_value_y_2;
2917  if (n_values > 1)
2918  test_value_y_2 = dim > 1 ?
2919  test_grads_value[3] * shapes[i1][0][1] :
2920  test_grads_value[3];
2921 
2922  Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
2923  Number2 *values_ptr_2 =
2924  n_values > 1 ?
2925  values_ptr + Utilities::fixed_power<dim>(n_shapes_runtime) :
2926  nullptr;
2927  for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
2928  {
2929  if (add)
2930  values_ptr[i0] += shapes[i0][0][0] * test_value_y;
2931  else
2932  values_ptr[i0] = shapes[i0][0][0] * test_value_y;
2933  values_ptr[i0] += shapes[i0][1][0] * test_grad_xy;
2934  if (n_values > 1)
2935  {
2936  if (add)
2937  values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2;
2938  else
2939  values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2;
2940  }
2941  }
2942  }
2943  i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
2944  }
2945  }
2946 
2947 
2948 
2953  template <int dim,
2954  typename Number,
2955  typename Number2,
2956  bool add,
2957  int n_values = 1>
2958  inline void
2960  const ::ndarray<Number, 2, dim> *shapes,
2961  const int n_shapes,
2962  const Number2 *value,
2963  const Tensor<1, dim, Number2> &gradient,
2964  Number2 *values)
2965  {
2966  static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
2967  static_assert(1 <= n_values && n_values <= 2,
2968  "Only n_values=1,2 implemented");
2969 
2970  // Note that 'add' is a template argument, so the compiler will remove
2971  // these checks
2972  if (dim == 0)
2973  {
2974  if (add)
2975  values[0] += value[0];
2976  else
2977  values[0] = value[0];
2978  if (n_values > 1)
2979  {
2980  if (add)
2981  values[1] += value[1];
2982  else
2983  values[1] = value[1];
2984  }
2985  return;
2986  }
2987 
2988  // Implement the transpose of the function above
2989  // as in evaluate, use `int` type to produce better code in this context
2990  std::array<Number2, 2 + n_values> test_grads_value;
2991  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2992  {
2993  // test grad x
2994  test_grads_value[0] =
2995  dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
2996  // test grad y
2997  test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] :
2998  (dim > 1 ? gradient[1] : Number2());
2999  // test value z
3000  test_grads_value[2] =
3001  dim > 2 ?
3002  (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
3003  value[0];
3004 
3005  if (n_values > 1)
3006  test_grads_value[3] =
3007  dim > 2 ? value[1] * shapes[i2][0][2] : value[1];
3008  // Generate separate code with known loop bounds for the most common
3009  // cases
3010  if (n_shapes == 2)
3011  do_apply_test_functions_xy<dim, 2, Number2, Number, add, n_values>(
3012  values, shapes, test_grads_value, n_shapes, i);
3013  else if (n_shapes == 3)
3014  do_apply_test_functions_xy<dim, 3, Number2, Number, add, n_values>(
3015  values, shapes, test_grads_value, n_shapes, i);
3016  else if (n_shapes == 4)
3017  do_apply_test_functions_xy<dim, 4, Number2, Number, add, n_values>(
3018  values, shapes, test_grads_value, n_shapes, i);
3019  else if (n_shapes == 5)
3020  do_apply_test_functions_xy<dim, 5, Number2, Number, add, n_values>(
3021  values, shapes, test_grads_value, n_shapes, i);
3022  else if (n_shapes == 6)
3023  do_apply_test_functions_xy<dim, 6, Number2, Number, add, n_values>(
3024  values, shapes, test_grads_value, n_shapes, i);
3025  else
3026  do_apply_test_functions_xy<dim, -1, Number2, Number, add, n_values>(
3027  values, shapes, test_grads_value, n_shapes, i);
3028  }
3029  }
3030 
3031 
3032 
3037  template <int dim,
3038  typename Number,
3039  typename Number2,
3040  bool add,
3041  int n_values = 1>
3042  inline void
3044  const Number2 *value,
3045  const Tensor<1, dim, Number2> &gradient,
3046  Number2 *values,
3047  const Point<dim, Number> &p)
3048  {
3049  static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
3050  static_assert(1 <= n_values && n_values <= 2,
3051  "Only n_values=1,2 implemented");
3052 
3053  // Note that 'add' is a template argument, so the compiler will remove
3054  // these checks
3055  if (dim == 0)
3056  {
3057  if (add)
3058  values[0] += value[0];
3059  else
3060  values[0] = value[0];
3061  if (n_values > 1)
3062  {
3063  if (add)
3064  values[1] += value[1];
3065  else
3066  values[1] = value[1];
3067  }
3068  }
3069  else if (dim == 1)
3070  {
3071  const Number2 difference = value[0] * p[0] + gradient[0];
3072  if (add)
3073  {
3074  values[0] += value[0] - difference;
3075  values[1] += difference;
3076  }
3077  else
3078  {
3079  values[0] = value[0] - difference;
3080  values[1] = difference;
3081  }
3082  if (n_values > 1)
3083  {
3084  const Number2 product = value[1] * p[0];
3085  if (add)
3086  {
3087  values[2] += value[1] - product;
3088  values[3] += product;
3089  }
3090  else
3091  {
3092  values[2] = value[1] - product;
3093  values[3] = product;
3094  }
3095  }
3096  }
3097  else if (dim == 2)
3098  {
3099  const Number2 test_value_y1 = value[0] * p[1] + gradient[1];
3100  const Number2 test_value_y0 = value[0] - test_value_y1;
3101  const Number2 test_grad_xy1 = gradient[0] * p[1];
3102  const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1;
3103  const Number2 value0 = p[0] * test_value_y0 + test_grad_xy0;
3104  const Number2 value1 = p[0] * test_value_y1 + test_grad_xy1;
3105 
3106  if (add)
3107  {
3108  values[0] += test_value_y0 - value0;
3109  values[1] += value0;
3110  values[2] += test_value_y1 - value1;
3111  values[3] += value1;
3112  }
3113  else
3114  {
3115  values[0] = test_value_y0 - value0;
3116  values[1] = value0;
3117  values[2] = test_value_y1 - value1;
3118  values[3] = value1;
3119  }
3120 
3121  if (n_values > 1)
3122  {
3123  const Number2 test_value_y1_2 = value[1] * p[1];
3124  const Number2 test_value_y0_2 = value[1] - test_value_y1_2;
3125  const Number2 value0_2 = p[0] * test_value_y0_2;
3126  const Number2 value1_2 = p[0] * test_value_y1_2;
3127 
3128  if (add)
3129  {
3130  values[4] += test_value_y0_2 - value0_2;
3131  values[5] += value0_2;
3132  values[6] += test_value_y1_2 - value1_2;
3133  values[7] += value1_2;
3134  }
3135  else
3136  {
3137  values[4] = test_value_y0_2 - value0_2;
3138  values[5] = value0_2;
3139  values[6] = test_value_y1_2 - value1_2;
3140  values[7] = value1_2;
3141  }
3142  }
3143  }
3144  else if (dim == 3)
3145  {
3146  Assert(n_values == 1, ExcNotImplemented());
3147 
3148  const Number2 test_value_z1 = value[0] * p[2] + gradient[2];
3149  const Number2 test_value_z0 = value[0] - test_value_z1;
3150  const Number2 test_grad_x1 = gradient[0] * p[2];
3151  const Number2 test_grad_x0 = gradient[0] - test_grad_x1;
3152  const Number2 test_grad_y1 = gradient[1] * p[2];
3153  const Number2 test_grad_y0 = gradient[1] - test_grad_y1;
3154 
3155  const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0;
3156  const Number2 test_value_y00 = test_value_z0 - test_value_y01;
3157  const Number2 test_grad_xy01 = test_grad_x0 * p[1];
3158  const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01;
3159  const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1;
3160  const Number2 test_value_y10 = test_value_z1 - test_value_y11;
3161  const Number2 test_grad_xy11 = test_grad_x1 * p[1];
3162  const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11;
3163 
3164  const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00;
3165  const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01;
3166  const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10;
3167  const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11;
3168 
3169  if (add)
3170  {
3171  values[0] += test_value_y00 - value00;
3172  values[1] += value00;
3173  values[2] += test_value_y01 - value01;
3174  values[3] += value01;
3175  values[4] += test_value_y10 - value10;
3176  values[5] += value10;
3177  values[6] += test_value_y11 - value11;
3178  values[7] += value11;
3179  }
3180  else
3181  {
3182  values[0] = test_value_y00 - value00;
3183  values[1] = value00;
3184  values[2] = test_value_y01 - value01;
3185  values[3] = value01;
3186  values[4] = test_value_y10 - value10;
3187  values[5] = value10;
3188  values[6] = test_value_y11 - value11;
3189  values[7] = value11;
3190  }
3191  }
3192  }
3193 
3194 
3195 
3201  template <bool is_linear,
3202  int dim,
3203  typename Number,
3204  typename Number2,
3205  int n_values = 1>
3206  inline void
3208  const ::ndarray<Number, 2, dim> *shapes,
3209  const unsigned int n_shapes,
3210  const Number2 *value,
3211  const Tensor<1, dim, Number2> &gradient,
3212  Number2 *values,
3213  const Point<dim, Number> &p,
3214  const bool do_add)
3215  {
3216  if (do_add)
3217  {
3218  if (is_linear)
3220  dim,
3221  Number,
3222  Number2,
3223  true,
3224  n_values>(value, gradient, values, p);
3225  else
3227  dim,
3228  Number,
3229  Number2,
3230  true,
3231  n_values>(shapes, n_shapes, value, gradient, values);
3232  }
3233  else
3234  {
3235  if (is_linear)
3237  dim,
3238  Number,
3239  Number2,
3240  false,
3241  n_values>(value, gradient, values, p);
3242  else
3244  dim,
3245  Number,
3246  Number2,
3247  false,
3248  n_values>(shapes, n_shapes, value, gradient, values);
3249  }
3250  }
3251 
3252 
3253 
3257  template <int dim, int length, typename Number2, typename Number, bool add>
3258  inline
3259 #ifndef DEBUG
3261 #endif
3262  void
3264  Number2 *values,
3265  const ::ndarray<Number, 2, dim> *shapes,
3266  const Number2 &test_value,
3267  const int n_shapes_runtime,
3268  int &i)
3269  {
3270  if (length > 0)
3271  {
3272  constexpr unsigned int array_size = length > 0 ? length : 1;
3273  std::array<Number, array_size> shape_values_x;
3274  for (unsigned int i1 = 0; i1 < array_size; ++i1)
3275  shape_values_x[i1] = shapes[i1][0][0];
3276  for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
3277  {
3278  const Number2 test_value_y =
3279  dim > 1 ? test_value * shapes[i1][0][1] : test_value;
3280 
3281  Number2 *values_ptr = values + i + i1 * length;
3282  for (unsigned int i0 = 0; i0 < length; ++i0)
3283  {
3284  if (add)
3285  values_ptr[i0] += shape_values_x[i0] * test_value_y;
3286  else
3287  values_ptr[i0] = shape_values_x[i0] * test_value_y;
3288  }
3289  }
3290  i += (dim > 1 ? length * length : length);
3291  }
3292  else
3293  {
3294  for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
3295  {
3296  const Number2 test_value_y =
3297  dim > 1 ? test_value * shapes[i1][0][1] : test_value;
3298 
3299  Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
3300  for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
3301  {
3302  if (add)
3303  values_ptr[i0] += shapes[i0][0][0] * test_value_y;
3304  else
3305  values_ptr[i0] = shapes[i0][0][0] * test_value_y;
3306  }
3307  }
3308  i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
3309  }
3310  }
3311 
3312 
3313 
3317  template <int dim, typename Number, typename Number2, bool add>
3318  inline void
3320  const ::ndarray<Number, 2, dim> *shapes,
3321  const int n_shapes,
3322  const Number2 &value,
3323  Number2 *values)
3324  {
3325  static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
3326 
3327  // as in evaluate, use `int` type to produce better code in this context
3328 
3329  if (dim == 0)
3330  {
3331  if (add)
3332  values[0] += value;
3333  else
3334  values[0] = value;
3335  return;
3336  }
3337 
3338  // Implement the transpose of the function above
3339  Number2 test_value;
3340  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3341  {
3342  // test value z
3343  test_value = dim > 2 ? value * shapes[i2][0][2] : value;
3344 
3345  // Generate separate code with known loop bounds for the most common
3346  // cases
3347  if (n_shapes == 2)
3348  do_apply_test_functions_xy_value<dim, 2, Number2, Number, add>(
3349  values, shapes, test_value, n_shapes, i);
3350  else if (n_shapes == 3)
3351  do_apply_test_functions_xy_value<dim, 3, Number2, Number, add>(
3352  values, shapes, test_value, n_shapes, i);
3353  else if (n_shapes == 4)
3354  do_apply_test_functions_xy_value<dim, 4, Number2, Number, add>(
3355  values, shapes, test_value, n_shapes, i);
3356  else if (n_shapes == 5)
3357  do_apply_test_functions_xy_value<dim, 5, Number2, Number, add>(
3358  values, shapes, test_value, n_shapes, i);
3359  else if (n_shapes == 6)
3360  do_apply_test_functions_xy_value<dim, 6, Number2, Number, add>(
3361  values, shapes, test_value, n_shapes, i);
3362  else
3363  do_apply_test_functions_xy_value<dim, -1, Number2, Number, add>(
3364  values, shapes, test_value, n_shapes, i);
3365  }
3366  }
3367 
3368 
3369 
3374  template <int dim, typename Number, typename Number2, bool add>
3375  inline void
3377  Number2 *values,
3378  const Point<dim, Number> &p)
3379  {
3380  static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
3381 
3382  if (dim == 0)
3383  {
3384  if (add)
3385  values[0] += value;
3386  else
3387  values[0] = value;
3388  }
3389  else if (dim == 1)
3390  {
3391  const auto x0 = 1. - p[0], x1 = p[0];
3392 
3393  if (add)
3394  {
3395  values[0] += value * x0;
3396  values[1] += value * x1;
3397  }
3398  else
3399  {
3400  values[0] = value * x0;
3401  values[1] = value * x1;
3402  }
3403  }
3404  else if (dim == 2)
3405  {
3406  const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
3407 
3408  const auto test_value_y0 = value * y0;
3409  const auto test_value_y1 = value * y1;
3410 
3411  if (add)
3412  {
3413  values[0] += x0 * test_value_y0;
3414  values[1] += x1 * test_value_y0;
3415  values[2] += x0 * test_value_y1;
3416  values[3] += x1 * test_value_y1;
3417  }
3418  else
3419  {
3420  values[0] = x0 * test_value_y0;
3421  values[1] = x1 * test_value_y0;
3422  values[2] = x0 * test_value_y1;
3423  values[3] = x1 * test_value_y1;
3424  }
3425  }
3426  else if (dim == 3)
3427  {
3428  const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
3429  z0 = 1. - p[2], z1 = p[2];
3430 
3431  const auto test_value_z0 = value * z0;
3432  const auto test_value_z1 = value * z1;
3433 
3434  const auto test_value_y00 = test_value_z0 * y0;
3435  const auto test_value_y01 = test_value_z0 * y1;
3436  const auto test_value_y10 = test_value_z1 * y0;
3437  const auto test_value_y11 = test_value_z1 * y1;
3438 
3439  if (add)
3440  {
3441  values[0] += x0 * test_value_y00;
3442  values[1] += x1 * test_value_y00;
3443  values[2] += x0 * test_value_y01;
3444  values[3] += x1 * test_value_y01;
3445  values[4] += x0 * test_value_y10;
3446  values[5] += x1 * test_value_y10;
3447  values[6] += x0 * test_value_y11;
3448  values[7] += x1 * test_value_y11;
3449  }
3450  else
3451  {
3452  values[0] = x0 * test_value_y00;
3453  values[1] = x1 * test_value_y00;
3454  values[2] = x0 * test_value_y01;
3455  values[3] = x1 * test_value_y01;
3456  values[4] = x0 * test_value_y10;
3457  values[5] = x1 * test_value_y10;
3458  values[6] = x0 * test_value_y11;
3459  values[7] = x1 * test_value_y11;
3460  }
3461  }
3462  }
3463 
3464 
3465 
3471  template <bool is_linear, int dim, typename Number, typename Number2>
3472  inline void
3473  integrate_tensor_product_value(const ::ndarray<Number, 2, dim> *shapes,
3474  const unsigned int n_shapes,
3475  const Number2 &value,
3476  Number2 *values,
3477  const Point<dim, Number> &p,
3478  const bool do_add)
3479  {
3480  if (do_add)
3481  {
3482  if (is_linear)
3484  Number,
3485  Number2,
3486  true>(value,
3487  values,
3488  p);
3489  else
3491  Number,
3492  Number2,
3493  true>(shapes,
3494  n_shapes,
3495  value,
3496  values);
3497  }
3498  else
3499  {
3500  if (is_linear)
3502  Number,
3503  Number2,
3504  false>(value,
3505  values,
3506  p);
3507  else
3509  Number,
3510  Number2,
3511  false>(shapes,
3512  n_shapes,
3513  value,
3514  values);
3515  }
3516  }
3517 
3518 
3519 
3520  template <int dim, int n_points_1d_template, typename Number>
3521  inline void
3522  weight_fe_q_dofs_by_entity(const Number *weights,
3523  const unsigned int n_components,
3524  const int n_points_1d_non_template,
3525  Number *data)
3526  {
3527  const int n_points_1d = n_points_1d_template != -1 ?
3528  n_points_1d_template :
3529  n_points_1d_non_template;
3530 
3531  Assert(n_points_1d > 0, ExcNotImplemented());
3532  Assert(n_points_1d < 100, ExcNotImplemented());
3533 
3534  unsigned int compressed_index[100];
3535  compressed_index[0] = 0;
3536  for (int i = 1; i < n_points_1d - 1; ++i)
3537  compressed_index[i] = 1;
3538  compressed_index[n_points_1d - 1] = 2;
3539 
3540  for (unsigned int c = 0; c < n_components; ++c)
3541  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3542  for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3543  {
3544  const unsigned int shift =
3545  9 * compressed_index[k] + 3 * compressed_index[j];
3546  data[0] *= weights[shift];
3547  // loop bound as int avoids compiler warnings in case n_points_1d
3548  // == 1 (polynomial degree 0)
3549  const Number weight = weights[shift + 1];
3550  for (int i = 1; i < n_points_1d - 1; ++i)
3551  data[i] *= weight;
3552  data[n_points_1d - 1] *= weights[shift + 2];
3553  data += n_points_1d;
3554  }
3555  }
3556 
3557 
3558  template <int dim, int n_points_1d_template, typename Number>
3559  inline void
3560  weight_fe_q_dofs_by_entity_shifted(const Number *weights,
3561  const unsigned int n_components,
3562  const int n_points_1d_non_template,
3563  Number *data)
3564  {
3565  const int n_points_1d = n_points_1d_template != -1 ?
3566  n_points_1d_template :
3567  n_points_1d_non_template;
3568 
3569  Assert((n_points_1d % 2) == 1,
3570  ExcMessage("The function can only with add number of points"));
3571  Assert(n_points_1d > 0, ExcNotImplemented());
3572  Assert(n_points_1d < 100, ExcNotImplemented());
3573 
3574  const unsigned int n_inside_1d = n_points_1d / 2;
3575 
3576  unsigned int compressed_index[100];
3577 
3578  unsigned int c = 0;
3579  for (int i = 0; i < n_inside_1d; ++i)
3580  compressed_index[c++] = 0;
3581  compressed_index[c++] = 1;
3582  for (int i = 0; i < n_inside_1d; ++i)
3583  compressed_index[c++] = 2;
3584 
3585  for (unsigned int c = 0; c < n_components; ++c)
3586  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3587  for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3588  {
3589  const unsigned int shift =
3590  9 * compressed_index[k] + 3 * compressed_index[j];
3591 
3592  unsigned int c = 0;
3593  const Number weight1 = weights[shift];
3594  for (int i = 0; i < n_inside_1d; ++i)
3595  data[c++] *= weight1;
3596  data[c++] *= weights[shift + 1];
3597  const Number weight2 = weights[shift + 2];
3598  for (int i = 0; i < n_inside_1d; ++i)
3599  data[c++] *= weight2;
3600  data += n_points_1d;
3601  }
3602  }
3603 
3604 
3605  template <int dim, int n_points_1d_template, typename Number>
3606  inline bool
3608  const unsigned int n_components,
3609  const int n_points_1d_non_template,
3610  Number *weights)
3611  {
3612  const int n_points_1d = n_points_1d_template != -1 ?
3613  n_points_1d_template :
3614  n_points_1d_non_template;
3615 
3616  Assert(n_points_1d > 0, ExcNotImplemented());
3617  Assert(n_points_1d < 100, ExcNotImplemented());
3618 
3619  unsigned int compressed_index[100];
3620  compressed_index[0] = 0;
3621  for (int i = 1; i < n_points_1d - 1; ++i)
3622  compressed_index[i] = 1;
3623  compressed_index[n_points_1d - 1] = 2;
3624 
3625  // Insert the number data into a storage position for weight,
3626  // ensuring that the array has either not been touched before
3627  // or the previous content is the same. In case the previous
3628  // content has a different value, we exit this function and
3629  // signal to outer functions that the compression was not possible.
3630  const auto check_and_set = [](Number &weight, const Number &data) {
3631  if (weight == Number(-1.0) || weight == data)
3632  {
3633  weight = data;
3634  return true; // success for the entry
3635  }
3636 
3637  return false; // failure for the entry
3638  };
3639 
3640  for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3641  weights[c] = Number(-1.0);
3642 
3643  for (unsigned int c = 0; c < n_components; ++c)
3644  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3645  for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
3646  ++j, data += n_points_1d)
3647  {
3648  const unsigned int shift =
3649  9 * compressed_index[k] + 3 * compressed_index[j];
3650 
3651  if (!check_and_set(weights[shift], data[0]))
3652  return false; // failure
3653 
3654  for (int i = 1; i < n_points_1d - 1; ++i)
3655  if (!check_and_set(weights[shift + 1], data[i]))
3656  return false; // failure
3657 
3658  if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
3659  return false; // failure
3660  }
3661 
3662  return true; // success
3663  }
3664 
3665 
3666  template <int dim, int n_points_1d_template, typename Number>
3667  inline bool
3669  const Number *data,
3670  const unsigned int n_components,
3671  const int n_points_1d_non_template,
3672  Number *weights)
3673  {
3674  const int n_points_1d = n_points_1d_template != -1 ?
3675  n_points_1d_template :
3676  n_points_1d_non_template;
3677 
3678  Assert((n_points_1d % 2) == 1,
3679  ExcMessage("The function can only with add number of points"));
3680  Assert(n_points_1d > 0, ExcNotImplemented());
3681  Assert(n_points_1d < 100, ExcNotImplemented());
3682 
3683  const unsigned int n_inside_1d = n_points_1d / 2;
3684 
3685  unsigned int compressed_index[100];
3686 
3687  unsigned int c = 0;
3688  for (int i = 0; i < n_inside_1d; ++i)
3689  compressed_index[c++] = 0;
3690  compressed_index[c++] = 1;
3691  for (int i = 0; i < n_inside_1d; ++i)
3692  compressed_index[c++] = 2;
3693 
3694  // Insert the number data into a storage position for weight,
3695  // ensuring that the array has either not been touched before
3696  // or the previous content is the same. In case the previous
3697  // content has a different value, we exit this function and
3698  // signal to outer functions that the compression was not possible.
3699  const auto check_and_set = [](Number &weight, const Number &data) {
3700  if (weight == Number(-1.0) || weight == data)
3701  {
3702  weight = data;
3703  return true; // success for the entry
3704  }
3705 
3706  return false; // failure for the entry
3707  };
3708 
3709  for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3710  weights[c] = Number(-1.0);
3711 
3712  for (unsigned int comp = 0; comp < n_components; ++comp)
3713  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3714  for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
3715  ++j, data += n_points_1d)
3716  {
3717  const unsigned int shift =
3718  9 * compressed_index[k] + 3 * compressed_index[j];
3719 
3720  unsigned int c = 0;
3721 
3722  for (int i = 0; i < n_inside_1d; ++i)
3723  if (!check_and_set(weights[shift], data[c++]))
3724  return false; // failure
3725 
3726  if (!check_and_set(weights[shift + 1], data[c++]))
3727  return false; // failure
3728 
3729  for (int i = 0; i < n_inside_1d; ++i)
3730  if (!check_and_set(weights[shift + 2], data[c++]))
3731  return false; // failure
3732  }
3733 
3734  return true; // success
3735  }
3736 
3737 
3738 } // end of namespace internal
3739 
3740 
3742 
3743 #endif
pointer data()
Definition: point.h:112
Definition: tensor.h:516
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_RESTRICT
Definition: config.h:111
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcMessage(std::string arg1)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:191
@ matrix
Contents is actually a matrix.
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:192
std::string to_string(const T &t)
Definition: patterns.h:2391
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
T fixed_power(const T t)
Definition: utilities.h:975
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, 2+n_values > do_interpolate_xy(const Number *values, const std::vector< unsigned int > &renumber, const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes_runtime, int &i)
void integrate_add_tensor_product_value_linear(const Number2 &value, Number2 *values, const Point< dim, Number > &p)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
ProductTypeNoPoint< Number, Number2 >::type do_interpolate_xy_value(const Number *values, const std::vector< unsigned int > &renumber, const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes_runtime, int &i)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool do_add)
void do_apply_test_functions_xy_value(Number2 *values, const ::ndarray< Number, 2, dim > *shapes, const Number2 &test_value, const int n_shapes_runtime, int &i)
void integrate_add_tensor_product_value_and_gradient_linear(const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p)
void integrate_add_tensor_product_value_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 &value, Number2 *values)
Tensor< 1, 1, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_higher_derivatives(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< 1, Number2 > &p, const std::vector< unsigned int > &renumber={})
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const Number *values, const Point< dim, Number2 > &p)
std::enable_if_t< contract_onto_face, void > interpolate_to_face(const Number2 *shape_values, const std::array< int, 2 > &n_blocks, const std::array< int, 2 > &steps, const Number *input, Number *DEAL_II_RESTRICT output, const int n_rows_runtime=0, const int stride_runtime=1)
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double >> &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
void integrate_add_tensor_product_value_and_gradient_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values)
void weight_fe_q_dofs_by_entity_shifted(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
void do_apply_test_functions_xy(Number2 *values, const ::ndarray< Number, 2, dim > *shapes, const std::array< Number2, 2+n_values > &test_grads_value, const int n_shapes_runtime, int &i)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
bool compute_weights_fe_q_dofs_by_entity(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
bool compute_weights_fe_q_dofs_by_entity_shifted(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool do_add)
static const unsigned int invalid_unsigned_int
Definition: types.h:221
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static void normal(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const bool add_into_result=false, const int subface_index_1d=0)
static void tangential(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const int subface_index_1d=0)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
void values(const Number in[], Number out[]) const
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
void values_one_line(const Number in[], Number out[]) const
void gradients(const Number in[], Number out[]) const
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int=0, const unsigned int=0)
static constexpr unsigned int n_rows_of_product
void hessians(const Number in[], Number out[]) const
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
static constexpr unsigned int n_columns_of_product
void gradients_one_line(const Number in[], Number out[]) const
void hessians_one_line(const Number in[], Number out[]) const
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition: shape_info.h:315
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
typename ProductType< Number, Number2 >::type type