deal.II version 9.8.0
\(\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\}}\)
Loading...
Searching...
No Matches
tensor_product_point_kernels.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_tensor_product_point_kernels_h
15#define dealii_matrix_free_tensor_product_point_kernels_h
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/ndarray.h>
20#include <deal.II/base/polynomial.h>
21#include <deal.II/base/utilities.h>
22
23#include <deal.II/matrix_free/shape_info.h>
24
25
26DEAL_II_NAMESPACE_OPEN
27
28
29
30namespace internal
31{
38 template <typename Number, typename Number2>
43
44 template <int dim, typename Number, typename Number2>
45 struct ProductTypeNoPoint<Point<dim, Number>, Number2>
46 {
47 using type = typename ProductType<Tensor<1, dim, Number>, Number2>::type;
48 };
49
50
51
56 template <int dim, typename Number>
57 inline void
58 compute_values_of_array(
59 ::ndarray<Number, 2, dim> *shapes,
60 const std::vector<Polynomials::Polynomial<double>> &poly,
61 const Point<dim, Number> &p,
62 const unsigned int derivative = 1)
63 {
64 const int n_shapes = poly.size();
65
66 // Evaluate 1d polynomials and their derivatives
67 std::array<Number, dim> point;
68 for (unsigned int d = 0; d < dim; ++d)
69 point[d] = p[d];
70 for (int i = 0; i < n_shapes; ++i)
71 poly[i].values_of_array(point, derivative, shapes[i].data());
72 }
73
74
75
79 template <typename Number>
80 inline void
81 compute_values_of_array(::ndarray<Number, 2, 0> *,
82 const std::vector<Polynomials::Polynomial<double>> &,
83 const Point<0, Number> &,
84 const unsigned int)
85 {
86 DEAL_II_ASSERT_UNREACHABLE();
87 }
88
89
90
98 template <int dim, typename Number>
99 inline void
100 compute_values_of_array_in_pairs(
101 ::ndarray<Number, 2, dim> *shapes,
102 const std::vector<Polynomials::Polynomial<double>> &poly,
103 const Point<dim, Number> &p0,
104 const Point<dim, Number> &p1)
105 {
106 // Use 'int' variable here to let the compiler apply additional
107 // optimizations, in particular regarding multiplications and additions in
108 // loop increments that are known not to overflow/wrap around (as is the
109 // case for unsigned int).
110 const int n_shapes = poly.size();
111
112 std::array<Number, 2 * dim> point, result;
113 for (unsigned int d = 0; d < dim; ++d)
114 point[d] = p0[d];
115 for (unsigned int d = 0; d < dim; ++d)
116 point[dim + d] = p1[d];
117 for (int i = 0; i < n_shapes; ++i)
118 {
119 poly[i].values_of_array(point, 0, &result);
120 for (unsigned int j = 0; j < 2; ++j)
121 for (unsigned int d = 0; d < dim; ++d)
122 shapes[j * n_shapes + i][0][d] = result[j * dim + d];
123 }
124 }
125
126
127
131 template <int dim,
132 int length,
133 typename Number2,
134 typename Number,
135 int n_values = 1,
136 bool do_renumber = true,
137 int stride = 1>
138 inline
139#ifndef DEBUG
140 DEAL_II_ALWAYS_INLINE
141#endif
142 std::array<typename ProductTypeNoPoint<Number, Number2>::type,
143 2 + n_values>
144 do_interpolate_xy(const Number *values,
145 const std::vector<unsigned int> &renumber,
146 const ::ndarray<Number2, 2, dim> *shapes,
147 const int n_shapes_runtime,
148 int &i)
149 {
150 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
151 static_assert(1 <= n_values && n_values <= 2,
152 "Only n_values=1,2 implemented");
153
154 const int n_shapes = length > 0 ? length : n_shapes_runtime;
155
156 // If n_values > 1, we want to interpolate from a second array,
157 // placed in the same array immediately after the main data. This
158 // is used to interpolate normal derivatives onto faces.
159 const Number *values_2 =
160 n_values > 1 ?
161 values + stride * (length > 0 ?
162 Utilities::pow(length, dim) :
163 Utilities::fixed_power<dim>(n_shapes_runtime)) :
164 nullptr;
165 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
166 std::array<Number3, 2 + n_values> result = {};
167 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
168 {
169 // Interpolation + derivative x direction
170 std::array<Number3, 1 + n_values> inner_result = {};
171
172 // Distinguish the inner loop based on whether we have a
173 // renumbering or not
174 if (do_renumber && !renumber.empty())
175 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
176 {
177 // gradient
178 inner_result[0] +=
179 shapes[i0][1][0] * values[renumber[i] * stride];
180 // values
181 inner_result[1] +=
182 shapes[i0][0][0] * values[renumber[i] * stride];
183 if (n_values > 1)
184 inner_result[2] +=
185 shapes[i0][0][0] * values_2[renumber[i] * stride];
186 }
187 else
188 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
189 {
190 // gradient
191 inner_result[0] += shapes[i0][1][0] * values[i * stride];
192 // values
193 inner_result[1] += shapes[i0][0][0] * values[i * stride];
194 if (n_values > 1)
195 inner_result[2] += shapes[i0][0][0] * values_2[i * stride];
196 }
197
198 if (dim > 1)
199 {
200 // Interpolation + derivative in y direction
201 // gradient
202 result[0] += inner_result[0] * shapes[i1][0][1];
203 result[1] += inner_result[1] * shapes[i1][1][1];
204 // values
205 result[2] += inner_result[1] * shapes[i1][0][1];
206 if (n_values > 1)
207 result[3] += inner_result[2] * shapes[i1][0][1];
208 }
209 else
210 {
211 // gradient
212 result[0] = inner_result[0];
213 // values
214 result[1] = inner_result[1];
215 if (n_values > 1)
216 result[2] = inner_result[2];
217 }
218 }
219 return result;
220 }
221
222
223
228 template <int dim,
229 typename Number,
230 typename Number2,
231 int n_values = 1,
232 bool do_renumber = true,
233 int stride = 1>
234 inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
235 dim + n_values>
236 evaluate_tensor_product_value_and_gradient_shapes(
237 const ::ndarray<Number2, 2, dim> *shapes,
238 const int n_shapes,
239 const Number *values,
240 const std::vector<unsigned int> &renumber = {})
241 {
242 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
243 static_assert(1 <= n_values && n_values <= 2,
244 "Only n_values=1,2 implemented");
245
246 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
247
248 std::array<Number3, dim + n_values> result = {};
249 if (dim == 0)
250 {
251 // We only need the interpolation of the value and normal derivatives on
252 // faces of a 1d element. As the interpolation is the value at the
253 // point, simply set the result vector accordingly.
254 result[0] = values[0];
255 if (n_values > 1)
256 result[1] = values[1 * stride];
257 return result;
258 }
259
260 // Go through the tensor product of shape functions and interpolate
261 // with optimal algorithm
262 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
263 {
264 std::array<Number3, 2 + n_values> inner_result;
265 // Generate separate code with known loop bounds for the most common
266 // cases
267 if (n_shapes == 2)
268 inner_result =
269 do_interpolate_xy<dim,
270 2,
271 Number2,
272 Number,
273 n_values,
274 do_renumber,
275 stride>(values, renumber, shapes, n_shapes, i);
276 else if (n_shapes == 3)
277 inner_result =
278 do_interpolate_xy<dim,
279 3,
280 Number2,
281 Number,
282 n_values,
283 do_renumber,
284 stride>(values, renumber, shapes, n_shapes, i);
285 else if (n_shapes == 4)
286 inner_result =
287 do_interpolate_xy<dim,
288 4,
289 Number2,
290 Number,
291 n_values,
292 do_renumber,
293 stride>(values, renumber, shapes, n_shapes, i);
294 else if (n_shapes == 5)
295 inner_result =
296 do_interpolate_xy<dim,
297 5,
298 Number2,
299 Number,
300 n_values,
301 do_renumber,
302 stride>(values, renumber, shapes, n_shapes, i);
303 else if (n_shapes == 6)
304 inner_result =
305 do_interpolate_xy<dim,
306 6,
307 Number2,
308 Number,
309 n_values,
310 do_renumber,
311 stride>(values, renumber, shapes, n_shapes, i);
312 else
313 inner_result =
314 do_interpolate_xy<dim,
315 -1,
316 Number2,
317 Number,
318 n_values,
319 do_renumber,
320 stride>(values, renumber, shapes, n_shapes, i);
321 if (dim == 3)
322 {
323 // derivative + interpolation in z direction
324 // gradient
325 result[0] += inner_result[0] * shapes[i2][0][2];
326 result[1] += inner_result[1] * shapes[i2][0][2];
327 result[2] += inner_result[2] * shapes[i2][1][2];
328 // values
329 result[3] += inner_result[2] * shapes[i2][0][2];
330 if (n_values > 1)
331 result[4] += inner_result[3] * shapes[i2][0][2];
332 }
333 else if (dim == 2)
334 {
335 // gradient
336 result[0] = inner_result[0];
337 result[1] = inner_result[1];
338 // values
339 result[2] = inner_result[2];
340 if (n_values > 1)
341 result[3] = inner_result[3];
342 }
343 else
344 {
345 // gradient
346 result[0] = inner_result[0];
347 // values
348 result[1] = inner_result[1];
349 if (n_values > 1)
350 result[2] = inner_result[2];
351 }
352 }
353
354 return result;
355 }
356
357
358
363 template <int dim,
364 typename Number,
365 typename Number2,
366 int n_values = 1,
367 int stride = 1>
368 inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
369 dim + n_values>
370 evaluate_tensor_product_value_and_gradient_linear(
371 const Number *values,
372 const Point<dim, Number2> &p)
373 {
374 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
375 static_assert(1 <= n_values && n_values <= 2,
376 "Only n_values=1,2 implemented");
377
378 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
379
380 // If n_values > 1, we want to interpolate from a second array,
381 // placed in the same array immediately after the main data. This
382 // is used to interpolate normal derivatives onto faces.
383
384 std::array<Number3, dim + n_values> result;
385 if (dim == 0)
386 {
387 // we only need the value on faces of a 1d element
388 result[0] = values[0];
389 if (n_values > 1)
390 result[1] = values[1 * stride];
391 }
392 else if (dim == 1)
393 {
394 // gradient
395 result[0] = Number3(values[stride] - values[0]);
396 // values
397 result[1] = Number3(values[0]) + p[0] * result[0];
398 if (n_values > 1)
399 result[2] = Number3(values[2 * stride]) +
400 p[0] * (values[3 * stride] - values[2 * stride]);
401 }
402 else if (dim == 2)
403 {
404 const Number3 val10 = Number3(values[stride] - values[0]);
405 const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
406 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
407 const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
408
409 // gradient
410 result[0] = val10 + p[1] * (val32 - val10);
411 result[1] = tmp1 - tmp0;
412
413 // values
414 result[2] = tmp0 + p[1] * result[1];
415
416 if (n_values > 1)
417 {
418 const Number3 tmp0_2 =
419 Number3(values[4 * stride]) +
420 p[0] * (values[5 * stride] - values[4 * stride]);
421 const Number3 tmp1_2 =
422 Number3(values[6 * stride]) +
423 p[0] * (values[7 * stride] - values[6 * stride]);
424 result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
425 }
426 }
427 else if (dim == 3)
428 {
429 const Number3 val10 = Number3(values[stride] - values[0]);
430 const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
431 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
432 const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
433 const Number3 tmp10 = tmp1 - tmp0;
434 const Number3 tmpy0 = tmp0 + p[1] * tmp10;
435
436 const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
437 const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
438 const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54;
439 const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76;
440 const Number3 tmp32 = tmp3 - tmp2;
441 const Number3 tmpy1 = tmp2 + p[1] * tmp32;
442
443 // gradient
444 result[2] = tmpy1 - tmpy0;
445 result[1] = tmp10 + p[2] * (tmp32 - tmp10);
446 const Number3 tmpz0 = val10 + p[1] * (val32 - val10);
447 result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0);
448
449 // value
450 result[3] = tmpy0 + p[2] * result[2];
451 Assert(n_values == 1, ExcNotImplemented());
452 }
453
454 return result;
455 }
456
457
458
493 template <int dim, typename Number, typename Number2>
494 inline std::pair<
495 typename ProductTypeNoPoint<Number, Number2>::type,
496 Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
497 evaluate_tensor_product_value_and_gradient(
498 const std::vector<Polynomials::Polynomial<double>> &poly,
499 const ArrayView<const Number> &values,
500 const Point<dim, Number2> &p,
501 const bool d_linear = false,
502 const std::vector<unsigned int> &renumber = {})
503 {
504 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
505
506 std::array<Number3, dim + 1> result;
507 if (d_linear)
508 {
509 result =
510 evaluate_tensor_product_value_and_gradient_linear(values.data(), p);
511 }
512 else
513 {
514 AssertIndexRange(poly.size(), 200);
515 std::array<::ndarray<Number2, 2, dim>, 200> shapes;
516 compute_values_of_array(shapes.data(), poly, p);
517 result = evaluate_tensor_product_value_and_gradient_shapes<dim,
518 Number,
519 Number2>(
520 shapes.data(), poly.size(), values.data(), renumber);
521 }
522 return std::make_pair(result[dim],
523 Tensor<1, dim, Number3>(
524 ArrayView<Number3>(result.data(), dim)));
525 }
526
527
528
529 template <int dim,
530 int length,
531 typename Number2,
532 typename Number,
533 bool do_renumber = true,
534 int stride = 1>
535 inline
536#ifndef DEBUG
537 DEAL_II_ALWAYS_INLINE
538#endif
539 typename ProductTypeNoPoint<Number, Number2>::type
540 do_interpolate_xy_value(const Number *values,
541 const std::vector<unsigned int> &renumber,
542 const ::ndarray<Number2, 2, dim> *shapes,
543 const int n_shapes_runtime,
544 int &i)
545 {
546 const int n_shapes = length > 0 ? length : n_shapes_runtime;
547 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
548 Number3 result = {};
549 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
550 {
551 // Interpolation x direction
552 Number3 value = {};
553
554 // Distinguish the inner loop based on whether we have a
555 // renumbering or not
556 if (do_renumber && !renumber.empty())
557 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
558 value += shapes[i0][0][0] * values[renumber[i] * stride];
559 else
560 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
561 value += shapes[i0][0][0] * values[i * stride];
562
563 if (dim > 1)
564 result += value * shapes[i1][0][1];
565 else
566 result = value;
567 }
568 return result;
569 }
570
571
572
573 template <int dim,
574 typename Number,
575 typename Number2,
576 bool do_renumber = true,
577 int stride = 1>
578 inline typename ProductTypeNoPoint<Number, Number2>::type
579 evaluate_tensor_product_value_shapes(
580 const ::ndarray<Number2, 2, dim> *shapes,
581 const int n_shapes,
582 const Number *values,
583 const std::vector<unsigned int> &renumber = {})
584 {
585 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
586
587 // we only need the value on faces of a 1d element
588 if (dim == 0)
589 {
590 return values[0];
591 }
592
593 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
594
595 // Go through the tensor product of shape functions and interpolate
596 // with optimal algorithm
597 Number3 result = {};
598 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
599 {
600 Number3 inner_result;
601 // Generate separate code with known loop bounds for the most common
602 // cases
603 if (n_shapes == 2)
604 inner_result = do_interpolate_xy_value<dim,
605 2,
606 Number2,
607 Number,
608 do_renumber,
609 stride>(
610 values, renumber, shapes, n_shapes, i);
611 else if (n_shapes == 3)
612 inner_result = do_interpolate_xy_value<dim,
613 3,
614 Number2,
615 Number,
616 do_renumber,
617 stride>(
618 values, renumber, shapes, n_shapes, i);
619 else if (n_shapes == 4)
620 inner_result = do_interpolate_xy_value<dim,
621 4,
622 Number2,
623 Number,
624 do_renumber,
625 stride>(
626 values, renumber, shapes, n_shapes, i);
627 else if (n_shapes == 5)
628 inner_result = do_interpolate_xy_value<dim,
629 5,
630 Number2,
631 Number,
632 do_renumber,
633 stride>(
634 values, renumber, shapes, n_shapes, i);
635 else if (n_shapes == 6)
636 inner_result = do_interpolate_xy_value<dim,
637 6,
638 Number2,
639 Number,
640 do_renumber,
641 stride>(
642 values, renumber, shapes, n_shapes, i);
643 else
644 inner_result = do_interpolate_xy_value<dim,
645 -1,
646 Number2,
647 Number,
648 do_renumber,
649 stride>(
650 values, renumber, shapes, n_shapes, i);
651 if (dim == 3)
652 {
653 // Interpolation + derivative in z direction
654 result += inner_result * shapes[i2][0][2];
655 }
656 else
657 {
658 result = inner_result;
659 }
660 }
661
662 return result;
663 }
664
665
666
667 template <int dim, typename Number, typename Number2, int stride = 1>
668 inline typename ProductTypeNoPoint<Number, Number2>::type
669 evaluate_tensor_product_value_linear(const Number *values,
670 const Point<dim, Number2> &p)
671 {
672 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
673
674 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
675
676 if (dim == 0)
677 {
678 // we only need the value on faces of a 1d element
679 return values[0];
680 }
681 else if (dim == 1)
682 {
683 return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]);
684 }
685 else if (dim == 2)
686 {
687 const Number3 val10 = Number3(values[stride] - values[0]);
688 const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
689 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
690 const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
691 return tmp0 + p[1] * (tmp1 - tmp0);
692 }
693 else if (dim == 3)
694 {
695 const Number3 val10 = Number3(values[stride] - values[0]);
696 const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
697 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
698 const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32;
699 const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
700
701 const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
702 const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
703 const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54;
704 const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76;
705 const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
706
707 return tmpy0 + p[2] * (tmpy1 - tmpy0);
708 }
709
710 // work around a compile error: missing return statement
711 return Number3();
712 }
713
714
715
716 template <int dim, typename Number, typename Number2>
717 inline typename ProductTypeNoPoint<Number, Number2>::type
718 evaluate_tensor_product_value(
719 const std::vector<Polynomials::Polynomial<double>> &poly,
720 const ArrayView<const Number> &values,
721 const Point<dim, Number2> &p,
722 const bool d_linear = false,
723 const std::vector<unsigned int> &renumber = {})
724 {
725 typename ProductTypeNoPoint<Number, Number2>::type result;
726 if (d_linear)
727 {
728 result = evaluate_tensor_product_value_linear(values.data(), p);
729 }
730 else
731 {
732 AssertIndexRange(poly.size(), 200);
733 std::array<::ndarray<Number2, 2, dim>, 200> shapes;
734 const int n_shapes = poly.size();
735 std::array<Number2, dim> point;
736 for (unsigned int d = 0; d < dim; ++d)
737 point[d] = p[d];
738 for (int i = 0; i < n_shapes; ++i)
739 poly[i].values_of_array(point, 0, shapes[i].data());
740 result = evaluate_tensor_product_value_shapes<dim, Number, Number2>(
741 shapes.data(), n_shapes, values.data(), renumber);
742 }
743 return result;
744 }
745
746
747
752 template <int derivative_order, typename Number, typename Number2>
753 inline Tensor<1, 1, typename ProductTypeNoPoint<Number, Number2>::type>
754 evaluate_tensor_product_higher_derivatives(
755 const std::vector<Polynomials::Polynomial<double>> &poly,
756 const ArrayView<const Number> &values,
757 const Point<1, Number2> &p,
758 const std::vector<unsigned int> &renumber = {})
759 {
760 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
761
762 const int n_shapes = poly.size();
763 AssertDimension(n_shapes, values.size());
764 Assert(renumber.empty() || renumber.size() == values.size(),
765 ExcDimensionMismatch(renumber.size(), values.size()));
766
767 std::array<Number2, derivative_order + 1> shapes;
768 Tensor<1, 1, Number3> result;
769 if (renumber.empty())
770 for (int i = 0; i < n_shapes; ++i)
771 {
772 poly[i].value(p[0], derivative_order, shapes.data());
773 result[0] += shapes[derivative_order] * values[i];
774 }
775 else
776 for (int i = 0; i < n_shapes; ++i)
777 {
778 poly[i].value(p[0], derivative_order, shapes.data());
779 result[0] += shapes[derivative_order] * values[renumber[i]];
780 }
781 return result;
782 }
783
784
785
790 template <int derivative_order, typename Number, typename Number2>
791 inline Tensor<1,
792 derivative_order + 1,
793 typename ProductTypeNoPoint<Number, Number2>::type>
794 evaluate_tensor_product_higher_derivatives(
795 const std::vector<Polynomials::Polynomial<double>> &poly,
796 const ArrayView<const Number> &values,
797 const Point<2, Number2> &p,
798 const std::vector<unsigned int> &renumber = {})
799 {
800 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
801 constexpr int dim = 2;
802
803 const int n_shapes = poly.size();
804 AssertDimension(Utilities::pow(n_shapes, 2), values.size());
805 Assert(renumber.empty() || renumber.size() == values.size(),
806 ExcDimensionMismatch(renumber.size(), values.size()));
807
808 AssertIndexRange(n_shapes, 100);
809 ::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
810 // Evaluate 1d polynomials and their derivatives
811 std::array<Number2, dim> point;
812 for (unsigned int d = 0; d < dim; ++d)
813 point[d] = p[d];
814 for (int i = 0; i < n_shapes; ++i)
815 poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
816
817 Tensor<1, derivative_order + 1, Number3> result;
818 for (int i1 = 0, i = 0; i1 < n_shapes; ++i1)
819 {
820 Tensor<1, derivative_order + 1, Number3> result_x;
821 if (renumber.empty())
822 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
823 for (unsigned int d = 0; d <= derivative_order; ++d)
824 result_x[d] += shapes[i0][d][0] * values[i];
825 else
826 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
827 for (unsigned int d = 0; d <= derivative_order; ++d)
828 result_x[d] += shapes[i0][d][0] * values[renumber[i]];
829
830 for (unsigned int d = 0; d <= derivative_order; ++d)
831 result[d] += shapes[i1][d][1] * result_x[derivative_order - d];
832 }
833 return result;
834 }
835
836
837
842 template <int derivative_order, typename Number, typename Number2>
843 inline Tensor<1,
844 ((derivative_order + 1) * (derivative_order + 2)) / 2,
845 typename ProductTypeNoPoint<Number, Number2>::type>
846 evaluate_tensor_product_higher_derivatives(
847 const std::vector<Polynomials::Polynomial<double>> &poly,
848 const ArrayView<const Number> &values,
849 const Point<3, Number2> &p,
850 const std::vector<unsigned int> &renumber = {})
851 {
852 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
853 constexpr int dim = 3;
854 constexpr int n_derivatives =
855 ((derivative_order + 1) * (derivative_order + 2)) / 2;
856
857 const int n_shapes = poly.size();
858 AssertDimension(Utilities::pow(n_shapes, 3), values.size());
859 Assert(renumber.empty() || renumber.size() == values.size(),
860 ExcDimensionMismatch(renumber.size(), values.size()));
861
862 AssertIndexRange(n_shapes, 100);
863 ::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
864 // Evaluate 1d polynomials and their derivatives
865 std::array<Number2, dim> point;
866 for (unsigned int d = 0; d < dim; ++d)
867 point[d] = p[d];
868 for (int i = 0; i < n_shapes; ++i)
869 poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
870
871 Tensor<1, n_derivatives, Number3> result;
872 for (int i2 = 0, i = 0; i2 < n_shapes; ++i2)
873 {
874 Tensor<1, n_derivatives, Number3> result_xy;
875 for (int i1 = 0; i1 < n_shapes; ++i1)
876 {
877 // apply x derivatives
878 Tensor<1, derivative_order + 1, Number3> result_x;
879 if (renumber.empty())
880 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
881 for (unsigned int d = 0; d <= derivative_order; ++d)
882 result_x[d] += shapes[i0][d][0] * values[i];
883 else
884 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
885 for (unsigned int d = 0; d <= derivative_order; ++d)
886 result_x[d] += shapes[i0][d][0] * values[renumber[i]];
887
888 // multiply by y derivatives, sorting them in upper triangular
889 // matrix, starting with highest global derivative order,
890 // decreasing the combined order of xy derivatives by one in each
891 // row, to be combined with z derivatives in the next step
892 for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
893 for (unsigned int e = d; e <= derivative_order; ++e, ++c)
894 result_xy[c] +=
895 shapes[i1][e - d][1] * result_x[derivative_order - e];
896 }
897
898 // multiply by z derivatives, starting with highest x derivative
899 for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
900 for (unsigned int e = d; e <= derivative_order; ++e, ++c)
901 result[c] += shapes[i2][d][2] * result_xy[c];
902 }
903 return result;
904 }
905
906
907
908 template <int dim, typename Number, typename Number2>
909 SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
910 evaluate_tensor_product_hessian(
911 const std::vector<Polynomials::Polynomial<double>> &poly,
912 const ArrayView<const Number> &values,
913 const Point<dim, Number2> &p,
914 const std::vector<unsigned int> &renumber = {})
915 {
916 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
917
918 const auto hessian =
919 evaluate_tensor_product_higher_derivatives<2>(poly, values, p, renumber);
920
921 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
922 SymmetricTensor<2, dim, Number3> result;
923 if (dim == 1)
924 result[0][0] = hessian[0];
925 else if (dim >= 2)
926 {
927 // derivatives in Hessian are xx, xy, yy, xz, yz, zz, so must re-order
928 // them for 3D
929 for (unsigned int d = 0, c = 0; d < 2; ++d)
930 for (unsigned int e = d; e < 2; ++e, ++c)
931 result[d][e] = hessian[c];
932 if (dim == 3)
933 {
934 for (unsigned int d = 0; d < 2; ++d)
935 result[d][2] = hessian[3 + d];
936 result[2][2] = hessian[5];
937 }
938 }
939
940 return result;
941 }
942
943
944
948 template <int dim,
949 int length,
950 typename Number2,
951 typename Number,
952 bool add,
953 int n_values = 1>
954 inline
955#ifndef DEBUG
956 DEAL_II_ALWAYS_INLINE
957#endif
958 void
959 do_apply_test_functions_xy(
960 Number2 *values,
961 const ::ndarray<Number, 2, dim> *shapes,
962 const std::array<Number2, 2 + n_values> &test_grads_value,
963 const int n_shapes_runtime,
964 int &i)
965 {
966 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
967 static_assert(1 <= n_values && n_values <= 2,
968 "Only n_values=1,2 implemented");
969
970 // Note that 'add' is a template argument, so the compiler will remove
971 // these checks
972 if (length > 0)
973 {
974 constexpr unsigned int array_size = length > 0 ? length : 1;
975 std::array<Number, array_size> shape_values_x;
976 std::array<Number, array_size> shape_derivs_x;
977 for (unsigned int j = 0; j < array_size; ++j)
978 {
979 shape_values_x[j] = shapes[j][0][0];
980 shape_derivs_x[j] = shapes[j][1][0];
981 }
982 for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
983 {
984 const Number2 test_value_y =
985 dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
986 test_grads_value[1] * shapes[i1][1][1]) :
987 test_grads_value[2];
988 const Number2 test_grad_xy =
989 dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
990 test_grads_value[0];
991 Number2 test_value_y_2;
992 if (n_values > 1)
993 test_value_y_2 = dim > 1 ?
994 test_grads_value[3] * shapes[i1][0][1] :
995 test_grads_value[3];
996
997 Number2 *values_ptr = values + i + i1 * length;
998 Number2 *values_ptr_2 =
999 n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr;
1000 for (int i0 = 0; i0 < length; ++i0)
1001 {
1002 if (add)
1003 values_ptr[i0] += shape_values_x[i0] * test_value_y;
1004 else
1005 values_ptr[i0] = shape_values_x[i0] * test_value_y;
1006 values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy;
1007 if (n_values > 1)
1008 {
1009 if (add)
1010 values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2;
1011 else
1012 values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2;
1013 }
1014 }
1015 }
1016 i += (dim > 1 ? length * length : length);
1017 }
1018 else
1019 {
1020 for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
1021 {
1022 const Number2 test_value_y =
1023 dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
1024 test_grads_value[1] * shapes[i1][1][1]) :
1025 test_grads_value[2];
1026 const Number2 test_grad_xy =
1027 dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
1028 test_grads_value[0];
1029 Number2 test_value_y_2;
1030 if (n_values > 1)
1031 test_value_y_2 = dim > 1 ?
1032 test_grads_value[3] * shapes[i1][0][1] :
1033 test_grads_value[3];
1034
1035 Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
1036 Number2 *values_ptr_2 =
1037 n_values > 1 ?
1038 values_ptr + Utilities::fixed_power<dim>(n_shapes_runtime) :
1039 nullptr;
1040 for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
1041 {
1042 if (add)
1043 values_ptr[i0] += shapes[i0][0][0] * test_value_y;
1044 else
1045 values_ptr[i0] = shapes[i0][0][0] * test_value_y;
1046 values_ptr[i0] += shapes[i0][1][0] * test_grad_xy;
1047 if (n_values > 1)
1048 {
1049 if (add)
1050 values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2;
1051 else
1052 values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2;
1053 }
1054 }
1055 }
1056 i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
1057 }
1058 }
1059
1060
1061
1066 template <int dim,
1067 typename Number,
1068 typename Number2,
1069 bool add,
1070 int n_values = 1>
1071 inline void
1072 integrate_add_tensor_product_value_and_gradient_shapes(
1073 const ::ndarray<Number, 2, dim> *shapes,
1074 const int n_shapes,
1075 const Number2 *value,
1076 const Tensor<1, dim, Number2> &gradient,
1077 Number2 *values)
1078 {
1079 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
1080 static_assert(1 <= n_values && n_values <= 2,
1081 "Only n_values=1,2 implemented");
1082
1083 // Note that 'add' is a template argument, so the compiler will remove
1084 // these checks
1085 if (dim == 0)
1086 {
1087 if (add)
1088 values[0] += value[0];
1089 else
1090 values[0] = value[0];
1091 if (n_values > 1)
1092 {
1093 if (add)
1094 values[1] += value[1];
1095 else
1096 values[1] = value[1];
1097 }
1098 return;
1099 }
1100
1101 // Implement the transpose of the function above
1102 // as in evaluate, use `int` type to produce better code in this context
1103 std::array<Number2, 2 + n_values> test_grads_value;
1104 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
1105 {
1106 // test grad x
1107 test_grads_value[0] =
1108 dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
1109 // test grad y
1110 test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] :
1111 (dim > 1 ? gradient[1] : Number2());
1112 // test value z
1113 test_grads_value[2] =
1114 dim > 2 ?
1115 (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
1116 value[0];
1117
1118 if (n_values > 1)
1119 test_grads_value[3] =
1120 dim > 2 ? value[1] * shapes[i2][0][2] : value[1];
1121 // Generate separate code with known loop bounds for the most common
1122 // cases
1123 if (n_shapes == 2)
1124 do_apply_test_functions_xy<dim, 2, Number2, Number, add, n_values>(
1125 values, shapes, test_grads_value, n_shapes, i);
1126 else if (n_shapes == 3)
1127 do_apply_test_functions_xy<dim, 3, Number2, Number, add, n_values>(
1128 values, shapes, test_grads_value, n_shapes, i);
1129 else if (n_shapes == 4)
1130 do_apply_test_functions_xy<dim, 4, Number2, Number, add, n_values>(
1131 values, shapes, test_grads_value, n_shapes, i);
1132 else if (n_shapes == 5)
1133 do_apply_test_functions_xy<dim, 5, Number2, Number, add, n_values>(
1134 values, shapes, test_grads_value, n_shapes, i);
1135 else if (n_shapes == 6)
1136 do_apply_test_functions_xy<dim, 6, Number2, Number, add, n_values>(
1137 values, shapes, test_grads_value, n_shapes, i);
1138 else
1139 do_apply_test_functions_xy<dim, -1, Number2, Number, add, n_values>(
1140 values, shapes, test_grads_value, n_shapes, i);
1141 }
1142 }
1143
1144
1145
1150 template <int dim,
1151 typename Number,
1152 typename Number2,
1153 bool add,
1154 int n_values = 1>
1155 inline void
1156 integrate_add_tensor_product_value_and_gradient_linear(
1157 const Number2 *value,
1158 const Tensor<1, dim, Number2> &gradient,
1159 Number2 *values,
1160 const Point<dim, Number> &p)
1161 {
1162 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
1163 static_assert(1 <= n_values && n_values <= 2,
1164 "Only n_values=1,2 implemented");
1165
1166 // Note that 'add' is a template argument, so the compiler will remove
1167 // these checks
1168 if (dim == 0)
1169 {
1170 if (add)
1171 values[0] += value[0];
1172 else
1173 values[0] = value[0];
1174 if (n_values > 1)
1175 {
1176 if (add)
1177 values[1] += value[1];
1178 else
1179 values[1] = value[1];
1180 }
1181 }
1182 else if (dim == 1)
1183 {
1184 const Number2 difference = value[0] * p[0] + gradient[0];
1185 if (add)
1186 {
1187 values[0] += value[0] - difference;
1188 values[1] += difference;
1189 }
1190 else
1191 {
1192 values[0] = value[0] - difference;
1193 values[1] = difference;
1194 }
1195 if (n_values > 1)
1196 {
1197 const Number2 product = value[1] * p[0];
1198 if (add)
1199 {
1200 values[2] += value[1] - product;
1201 values[3] += product;
1202 }
1203 else
1204 {
1205 values[2] = value[1] - product;
1206 values[3] = product;
1207 }
1208 }
1209 }
1210 else if (dim == 2)
1211 {
1212 const Number2 test_value_y1 = value[0] * p[1] + gradient[1];
1213 const Number2 test_value_y0 = value[0] - test_value_y1;
1214 const Number2 test_grad_xy1 = gradient[0] * p[1];
1215 const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1;
1216 const Number2 value0 = p[0] * test_value_y0 + test_grad_xy0;
1217 const Number2 value1 = p[0] * test_value_y1 + test_grad_xy1;
1218
1219 if (add)
1220 {
1221 values[0] += test_value_y0 - value0;
1222 values[1] += value0;
1223 values[2] += test_value_y1 - value1;
1224 values[3] += value1;
1225 }
1226 else
1227 {
1228 values[0] = test_value_y0 - value0;
1229 values[1] = value0;
1230 values[2] = test_value_y1 - value1;
1231 values[3] = value1;
1232 }
1233
1234 if (n_values > 1)
1235 {
1236 const Number2 test_value_y1_2 = value[1] * p[1];
1237 const Number2 test_value_y0_2 = value[1] - test_value_y1_2;
1238 const Number2 value0_2 = p[0] * test_value_y0_2;
1239 const Number2 value1_2 = p[0] * test_value_y1_2;
1240
1241 if (add)
1242 {
1243 values[4] += test_value_y0_2 - value0_2;
1244 values[5] += value0_2;
1245 values[6] += test_value_y1_2 - value1_2;
1246 values[7] += value1_2;
1247 }
1248 else
1249 {
1250 values[4] = test_value_y0_2 - value0_2;
1251 values[5] = value0_2;
1252 values[6] = test_value_y1_2 - value1_2;
1253 values[7] = value1_2;
1254 }
1255 }
1256 }
1257 else if (dim == 3)
1258 {
1259 Assert(n_values == 1, ExcNotImplemented());
1260
1261 const Number2 test_value_z1 = value[0] * p[2] + gradient[2];
1262 const Number2 test_value_z0 = value[0] - test_value_z1;
1263 const Number2 test_grad_x1 = gradient[0] * p[2];
1264 const Number2 test_grad_x0 = gradient[0] - test_grad_x1;
1265 const Number2 test_grad_y1 = gradient[1] * p[2];
1266 const Number2 test_grad_y0 = gradient[1] - test_grad_y1;
1267
1268 const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0;
1269 const Number2 test_value_y00 = test_value_z0 - test_value_y01;
1270 const Number2 test_grad_xy01 = test_grad_x0 * p[1];
1271 const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01;
1272 const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1;
1273 const Number2 test_value_y10 = test_value_z1 - test_value_y11;
1274 const Number2 test_grad_xy11 = test_grad_x1 * p[1];
1275 const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11;
1276
1277 const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00;
1278 const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01;
1279 const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10;
1280 const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11;
1281
1282 if (add)
1283 {
1284 values[0] += test_value_y00 - value00;
1285 values[1] += value00;
1286 values[2] += test_value_y01 - value01;
1287 values[3] += value01;
1288 values[4] += test_value_y10 - value10;
1289 values[5] += value10;
1290 values[6] += test_value_y11 - value11;
1291 values[7] += value11;
1292 }
1293 else
1294 {
1295 values[0] = test_value_y00 - value00;
1296 values[1] = value00;
1297 values[2] = test_value_y01 - value01;
1298 values[3] = value01;
1299 values[4] = test_value_y10 - value10;
1300 values[5] = value10;
1301 values[6] = test_value_y11 - value11;
1302 values[7] = value11;
1303 }
1304 }
1305 }
1306
1307
1308
1314 template <bool is_linear,
1315 int dim,
1316 typename Number,
1317 typename Number2,
1318 int n_values = 1>
1319 inline void
1320 integrate_tensor_product_value_and_gradient(
1321 const ::ndarray<Number, 2, dim> *shapes,
1322 const unsigned int n_shapes,
1323 const Number2 *value,
1324 const Tensor<1, dim, Number2> &gradient,
1325 Number2 *values,
1326 const Point<dim, Number> &p,
1327 const bool do_add)
1328 {
1329 if (do_add)
1330 {
1331 if (is_linear)
1332 internal::integrate_add_tensor_product_value_and_gradient_linear<
1333 dim,
1334 Number,
1335 Number2,
1336 true,
1337 n_values>(value, gradient, values, p);
1338 else
1339 internal::integrate_add_tensor_product_value_and_gradient_shapes<
1340 dim,
1341 Number,
1342 Number2,
1343 true,
1344 n_values>(shapes, n_shapes, value, gradient, values);
1345 }
1346 else
1347 {
1348 if (is_linear)
1349 internal::integrate_add_tensor_product_value_and_gradient_linear<
1350 dim,
1351 Number,
1352 Number2,
1353 false,
1354 n_values>(value, gradient, values, p);
1355 else
1356 internal::integrate_add_tensor_product_value_and_gradient_shapes<
1357 dim,
1358 Number,
1359 Number2,
1360 false,
1361 n_values>(shapes, n_shapes, value, gradient, values);
1362 }
1363 }
1364
1365
1366
1370 template <int dim, int length, typename Number2, typename Number, bool add>
1371 inline
1372#ifndef DEBUG
1373 DEAL_II_ALWAYS_INLINE
1374#endif
1375 void
1376 do_apply_test_functions_xy_value(
1377 Number2 *values,
1378 const ::ndarray<Number, 2, dim> *shapes,
1379 const Number2 &test_value,
1380 const int n_shapes_runtime,
1381 int &i)
1382 {
1383 if (length > 0)
1384 {
1385 constexpr unsigned int array_size = length > 0 ? length : 1;
1386 std::array<Number, array_size> shape_values_x;
1387 for (unsigned int i1 = 0; i1 < array_size; ++i1)
1388 shape_values_x[i1] = shapes[i1][0][0];
1389 for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
1390 {
1391 const Number2 test_value_y =
1392 dim > 1 ? test_value * shapes[i1][0][1] : test_value;
1393
1394 Number2 *values_ptr = values + i + i1 * length;
1395 for (int i0 = 0; i0 < length; ++i0)
1396 {
1397 if (add)
1398 values_ptr[i0] += shape_values_x[i0] * test_value_y;
1399 else
1400 values_ptr[i0] = shape_values_x[i0] * test_value_y;
1401 }
1402 }
1403 i += (dim > 1 ? length * length : length);
1404 }
1405 else
1406 {
1407 for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
1408 {
1409 const Number2 test_value_y =
1410 dim > 1 ? test_value * shapes[i1][0][1] : test_value;
1411
1412 Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
1413 for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
1414 {
1415 if (add)
1416 values_ptr[i0] += shapes[i0][0][0] * test_value_y;
1417 else
1418 values_ptr[i0] = shapes[i0][0][0] * test_value_y;
1419 }
1420 }
1421 i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
1422 }
1423 }
1424
1425
1426
1430 template <int dim, typename Number, typename Number2, bool add>
1431 inline void
1432 integrate_add_tensor_product_value_shapes(
1433 const ::ndarray<Number, 2, dim> *shapes,
1434 const int n_shapes,
1435 const Number2 &value,
1436 Number2 *values)
1437 {
1438 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
1439
1440 // as in evaluate, use `int` type to produce better code in this context
1441
1442 if (dim == 0)
1443 {
1444 if (add)
1445 values[0] += value;
1446 else
1447 values[0] = value;
1448 return;
1449 }
1450
1451 // Implement the transpose of the function above
1452 Number2 test_value;
1453 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
1454 {
1455 // test value z
1456 test_value = dim > 2 ? value * shapes[i2][0][2] : value;
1457
1458 // Generate separate code with known loop bounds for the most common
1459 // cases
1460 if (n_shapes == 2)
1461 do_apply_test_functions_xy_value<dim, 2, Number2, Number, add>(
1462 values, shapes, test_value, n_shapes, i);
1463 else if (n_shapes == 3)
1464 do_apply_test_functions_xy_value<dim, 3, Number2, Number, add>(
1465 values, shapes, test_value, n_shapes, i);
1466 else if (n_shapes == 4)
1467 do_apply_test_functions_xy_value<dim, 4, Number2, Number, add>(
1468 values, shapes, test_value, n_shapes, i);
1469 else if (n_shapes == 5)
1470 do_apply_test_functions_xy_value<dim, 5, Number2, Number, add>(
1471 values, shapes, test_value, n_shapes, i);
1472 else if (n_shapes == 6)
1473 do_apply_test_functions_xy_value<dim, 6, Number2, Number, add>(
1474 values, shapes, test_value, n_shapes, i);
1475 else
1476 do_apply_test_functions_xy_value<dim, -1, Number2, Number, add>(
1477 values, shapes, test_value, n_shapes, i);
1478 }
1479 }
1480
1481
1482
1487 template <int dim, typename Number, typename Number2, bool add>
1488 inline void
1489 integrate_add_tensor_product_value_linear(const Number2 &value,
1490 Number2 *values,
1491 const Point<dim, Number> &p)
1492 {
1493 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
1494
1495 if (dim == 0)
1496 {
1497 if (add)
1498 values[0] += value;
1499 else
1500 values[0] = value;
1501 }
1502 else if (dim == 1)
1503 {
1504 const auto x0 = Number(1.) - p[0], x1 = p[0];
1505
1506 if (add)
1507 {
1508 values[0] += value * x0;
1509 values[1] += value * x1;
1510 }
1511 else
1512 {
1513 values[0] = value * x0;
1514 values[1] = value * x1;
1515 }
1516 }
1517 else if (dim == 2)
1518 {
1519 const auto x0 = Number(1.) - p[0], x1 = p[0], y0 = Number(1.) - p[1],
1520 y1 = p[1];
1521
1522 const auto test_value_y0 = value * y0;
1523 const auto test_value_y1 = value * y1;
1524
1525 if (add)
1526 {
1527 values[0] += x0 * test_value_y0;
1528 values[1] += x1 * test_value_y0;
1529 values[2] += x0 * test_value_y1;
1530 values[3] += x1 * test_value_y1;
1531 }
1532 else
1533 {
1534 values[0] = x0 * test_value_y0;
1535 values[1] = x1 * test_value_y0;
1536 values[2] = x0 * test_value_y1;
1537 values[3] = x1 * test_value_y1;
1538 }
1539 }
1540 else if (dim == 3)
1541 {
1542 const auto x0 = Number(1.) - p[0], x1 = p[0], y0 = Number(1.) - p[1],
1543 y1 = p[1], z0 = Number(1.) - p[2], z1 = p[2];
1544
1545 const auto test_value_z0 = value * z0;
1546 const auto test_value_z1 = value * z1;
1547
1548 const auto test_value_y00 = test_value_z0 * y0;
1549 const auto test_value_y01 = test_value_z0 * y1;
1550 const auto test_value_y10 = test_value_z1 * y0;
1551 const auto test_value_y11 = test_value_z1 * y1;
1552
1553 if (add)
1554 {
1555 values[0] += x0 * test_value_y00;
1556 values[1] += x1 * test_value_y00;
1557 values[2] += x0 * test_value_y01;
1558 values[3] += x1 * test_value_y01;
1559 values[4] += x0 * test_value_y10;
1560 values[5] += x1 * test_value_y10;
1561 values[6] += x0 * test_value_y11;
1562 values[7] += x1 * test_value_y11;
1563 }
1564 else
1565 {
1566 values[0] = x0 * test_value_y00;
1567 values[1] = x1 * test_value_y00;
1568 values[2] = x0 * test_value_y01;
1569 values[3] = x1 * test_value_y01;
1570 values[4] = x0 * test_value_y10;
1571 values[5] = x1 * test_value_y10;
1572 values[6] = x0 * test_value_y11;
1573 values[7] = x1 * test_value_y11;
1574 }
1575 }
1576 }
1577
1578
1579
1585 template <bool is_linear, int dim, typename Number, typename Number2>
1586 inline void
1587 integrate_tensor_product_value(const ::ndarray<Number, 2, dim> *shapes,
1588 const unsigned int n_shapes,
1589 const Number2 &value,
1590 Number2 *values,
1591 const Point<dim, Number> &p,
1592 const bool do_add)
1593 {
1594 if (do_add)
1595 {
1596 if (is_linear)
1597 internal::integrate_add_tensor_product_value_linear<dim,
1598 Number,
1599 Number2,
1600 true>(value,
1601 values,
1602 p);
1603 else
1604 internal::integrate_add_tensor_product_value_shapes<dim,
1605 Number,
1606 Number2,
1607 true>(shapes,
1608 n_shapes,
1609 value,
1610 values);
1611 }
1612 else
1613 {
1614 if (is_linear)
1615 internal::integrate_add_tensor_product_value_linear<dim,
1616 Number,
1617 Number2,
1618 false>(value,
1619 values,
1620 p);
1621 else
1622 internal::integrate_add_tensor_product_value_shapes<dim,
1623 Number,
1624 Number2,
1625 false>(shapes,
1626 n_shapes,
1627 value,
1628 values);
1629 }
1630 }
1631} // end of namespace internal
1632
1633
1634DEAL_II_NAMESPACE_CLOSE
1635
1636#endif
Point
Definition point.h:111
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:165
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exception_macros.h:781
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
value
const bool IsBlockVector< VectorType >::value
Definition block_vector_base.h:91
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:210
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
PlasticityLab::Number
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
Definition ALE_Finite_Strain_Plasticity.h:8749
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::evaluate_tensor_product_value_and_gradient_shapes
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={})
Definition tensor_product_point_kernels.h:236
internal::do_interpolate_xy
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)
Definition tensor_product_point_kernels.h:144
internal::integrate_add_tensor_product_value_linear
void integrate_add_tensor_product_value_linear(const Number2 &value, Number2 *values, const Point< dim, Number > &p)
Definition tensor_product_point_kernels.h:1489
internal::integrate_tensor_product_value_and_gradient
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)
Definition tensor_product_point_kernels.h:1320
internal::do_apply_test_functions_xy_value
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)
Definition tensor_product_point_kernels.h:1376
internal::compute_values_of_array_in_pairs
void compute_values_of_array_in_pairs(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p0, const Point< dim, Number > &p1)
Definition tensor_product_point_kernels.h:100
internal::integrate_add_tensor_product_value_and_gradient_linear
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)
Definition tensor_product_point_kernels.h:1156
internal::integrate_add_tensor_product_value_shapes
void integrate_add_tensor_product_value_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 &value, Number2 *values)
Definition tensor_product_point_kernels.h:1432
internal::evaluate_tensor_product_hessian
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
Definition tensor_product_point_kernels.h:910
internal::evaluate_tensor_product_value_shapes
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={})
Definition tensor_product_point_kernels.h:579
internal::compute_values_of_array
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)
Definition tensor_product_point_kernels.h:58
internal::integrate_add_tensor_product_value_and_gradient_shapes
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)
Definition tensor_product_point_kernels.h:1072
internal::evaluate_tensor_product_value_and_gradient
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 ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
Definition tensor_product_point_kernels.h:497
internal::do_apply_test_functions_xy
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)
Definition tensor_product_point_kernels.h:959
internal::evaluate_tensor_product_value_linear
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
Definition tensor_product_point_kernels.h:669
internal::evaluate_tensor_product_value_and_gradient_linear
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)
Definition tensor_product_point_kernels.h:370
internal::evaluate_tensor_product_value
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
Definition tensor_product_point_kernels.h:718
internal::evaluate_tensor_product_higher_derivatives
Tensor< 1, 1, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_higher_derivatives(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< 1, Number2 > &p, const std::vector< unsigned int > &renumber={})
Definition tensor_product_point_kernels.h:754
internal::integrate_tensor_product_value
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)
Definition tensor_product_point_kernels.h:1587
internal::do_interpolate_xy_value
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)
Definition tensor_product_point_kernels.h:540
ndarray
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:105
shape_info.h
ProductType::type
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
Definition template_constraints.h:451
internal::ProductTypeNoPoint< Point< dim, Number >, Number2 >::type
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
Definition tensor_product_point_kernels.h:47
internal::ProductTypeNoPoint::type
typename ProductType< Number, Number2 >::type type
Definition tensor_product_point_kernels.h:41
Generated by doxygen 1.16.1