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