Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+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
function_lib.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 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
18#include <deal.II/base/point.h>
20#include <deal.II/base/tensor.h>
21
22#include <deal.II/lac/vector.h>
23
24#include <cmath>
25
27
28
29namespace Functions
30{
31 template <int dim>
32 double
33 SquareFunction<dim>::value(const Point<dim> &p, const unsigned int) const
34 {
35 return p.square();
36 }
37
38
39 template <int dim>
40 void
42 Vector<double> &values) const
43 {
44 AssertDimension(values.size(), 1);
45 values(0) = p.square();
46 }
47
48
49 template <int dim>
50 void
51 SquareFunction<dim>::value_list(const std::vector<Point<dim>> &points,
52 std::vector<double> &values,
53 const unsigned int) const
54 {
55 Assert(values.size() == points.size(),
56 ExcDimensionMismatch(values.size(), points.size()));
57
58 for (unsigned int i = 0; i < points.size(); ++i)
59 {
60 const Point<dim> &p = points[i];
61 values[i] = p.square();
62 }
63 }
64
65
66 template <int dim>
67 double
68 SquareFunction<dim>::laplacian(const Point<dim> &, const unsigned int) const
69 {
70 return 2 * dim;
71 }
72
73
74 template <int dim>
75 void
77 std::vector<double> &values,
78 const unsigned int) const
79 {
80 Assert(values.size() == points.size(),
81 ExcDimensionMismatch(values.size(), points.size()));
82
83 for (unsigned int i = 0; i < points.size(); ++i)
84 values[i] = 2 * dim;
85 }
86
87
88
89 template <int dim>
91 SquareFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
92 {
93 return p * 2.;
94 }
95
96
97 template <int dim>
98 void
100 const Point<dim> &p,
101 std::vector<Tensor<1, dim>> &values) const
102 {
103 AssertDimension(values.size(), 1);
104 values[0] = p * 2.;
105 }
106
107
108
109 template <int dim>
110 void
112 std::vector<Tensor<1, dim>> &gradients,
113 const unsigned int) const
114 {
115 Assert(gradients.size() == points.size(),
116 ExcDimensionMismatch(gradients.size(), points.size()));
117
118 for (unsigned int i = 0; i < points.size(); ++i)
119 gradients[i] = static_cast<Tensor<1, dim>>(points[i]) * 2;
120 }
121
122
123 //--------------------------------------------------------------------
124
125
126 template <int dim>
127 double
128 Q1WedgeFunction<dim>::value(const Point<dim> &p, const unsigned int) const
129 {
130 Assert(dim >= 2, ExcInternalError());
131 return p[0] * p[1];
132 }
133
134
135
136 template <int dim>
137 void
139 std::vector<double> &values,
140 const unsigned int) const
141 {
142 Assert(dim >= 2, ExcInternalError());
143 Assert(values.size() == points.size(),
144 ExcDimensionMismatch(values.size(), points.size()));
145
146 for (unsigned int i = 0; i < points.size(); ++i)
147 {
148 const Point<dim> &p = points[i];
149 values[i] = p[0] * p[1];
150 }
151 }
152
153
154 template <int dim>
155 void
157 const std::vector<Point<dim>> &points,
158 std::vector<Vector<double>> &values) const
159 {
160 Assert(dim >= 2, ExcInternalError());
161 Assert(values.size() == points.size(),
162 ExcDimensionMismatch(values.size(), points.size()));
163 Assert(values[0].size() == 1, ExcDimensionMismatch(values[0].size(), 1));
164
165 for (unsigned int i = 0; i < points.size(); ++i)
166 {
167 const Point<dim> &p = points[i];
168 values[i](0) = p[0] * p[1];
169 }
170 }
171
172
173 template <int dim>
174 double
175 Q1WedgeFunction<dim>::laplacian(const Point<dim> &, const unsigned int) const
176 {
177 Assert(dim >= 2, ExcInternalError());
178 return 0.;
179 }
180
181
182 template <int dim>
183 void
185 std::vector<double> &values,
186 const unsigned int) const
187 {
188 Assert(dim >= 2, ExcInternalError());
189 Assert(values.size() == points.size(),
190 ExcDimensionMismatch(values.size(), points.size()));
191
192 for (unsigned int i = 0; i < points.size(); ++i)
193 values[i] = 0.;
194 }
195
196
197
198 template <int dim>
200 Q1WedgeFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
201 {
202 Assert(dim >= 2, ExcInternalError());
203 Tensor<1, dim> erg;
204 erg[0] = p[1];
205 erg[1] = p[0];
206 return erg;
207 }
208
209
210
211 template <int dim>
212 void
214 std::vector<Tensor<1, dim>> &gradients,
215 const unsigned int) const
216 {
217 Assert(dim >= 2, ExcInternalError());
218 Assert(gradients.size() == points.size(),
219 ExcDimensionMismatch(gradients.size(), points.size()));
220
221 for (unsigned int i = 0; i < points.size(); ++i)
222 {
223 gradients[i][0] = points[i][1];
224 gradients[i][1] = points[i][0];
225 }
226 }
227
228
229 template <int dim>
230 void
232 const std::vector<Point<dim>> &points,
233 std::vector<std::vector<Tensor<1, dim>>> &gradients) const
234 {
235 Assert(dim >= 2, ExcInternalError());
236 Assert(gradients.size() == points.size(),
237 ExcDimensionMismatch(gradients.size(), points.size()));
238 Assert(gradients[0].size() == 1,
239 ExcDimensionMismatch(gradients[0].size(), 1));
240
241 for (unsigned int i = 0; i < points.size(); ++i)
242 {
243 gradients[i][0][0] = points[i][1];
244 gradients[i][0][1] = points[i][0];
245 }
246 }
247
248
249 //--------------------------------------------------------------------
250
251
252 template <int dim>
254 : offset(offset)
255 {}
256
257
258 template <int dim>
259 double
260 PillowFunction<dim>::value(const Point<dim> &p, const unsigned int) const
261 {
262 switch (dim)
263 {
264 case 1:
265 return 1. - p[0] * p[0] + offset;
266 case 2:
267 return (1. - p[0] * p[0]) * (1. - p[1] * p[1]) + offset;
268 case 3:
269 return (1. - p[0] * p[0]) * (1. - p[1] * p[1]) * (1. - p[2] * p[2]) +
270 offset;
271 default:
273 }
274 return 0.;
275 }
276
277 template <int dim>
278 void
279 PillowFunction<dim>::value_list(const std::vector<Point<dim>> &points,
280 std::vector<double> &values,
281 const unsigned int) const
282 {
283 Assert(values.size() == points.size(),
284 ExcDimensionMismatch(values.size(), points.size()));
285
286 for (unsigned int i = 0; i < points.size(); ++i)
287 {
288 const Point<dim> &p = points[i];
289 switch (dim)
290 {
291 case 1:
292 values[i] = 1. - p[0] * p[0] + offset;
293 break;
294 case 2:
295 values[i] = (1. - p[0] * p[0]) * (1. - p[1] * p[1]) + offset;
296 break;
297 case 3:
298 values[i] =
299 (1. - p[0] * p[0]) * (1. - p[1] * p[1]) * (1. - p[2] * p[2]) +
300 offset;
301 break;
302 default:
304 }
305 }
306 }
307
308
309
310 template <int dim>
311 double
312 PillowFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
313 {
314 switch (dim)
315 {
316 case 1:
317 return -2.;
318 case 2:
319 return -2. * ((1. - p[0] * p[0]) + (1. - p[1] * p[1]));
320 case 3:
321 return -2. * ((1. - p[0] * p[0]) * (1. - p[1] * p[1]) +
322 (1. - p[1] * p[1]) * (1. - p[2] * p[2]) +
323 (1. - p[2] * p[2]) * (1. - p[0] * p[0]));
324 default:
326 }
327 return 0.;
328 }
329
330 template <int dim>
331 void
333 std::vector<double> &values,
334 const unsigned int) const
335 {
336 Assert(values.size() == points.size(),
337 ExcDimensionMismatch(values.size(), points.size()));
338
339 for (unsigned int i = 0; i < points.size(); ++i)
340 {
341 const Point<dim> &p = points[i];
342 switch (dim)
343 {
344 case 1:
345 values[i] = -2.;
346 break;
347 case 2:
348 values[i] = -2. * ((1. - p[0] * p[0]) + (1. - p[1] * p[1]));
349 break;
350 case 3:
351 values[i] = -2. * ((1. - p[0] * p[0]) * (1. - p[1] * p[1]) +
352 (1. - p[1] * p[1]) * (1. - p[2] * p[2]) +
353 (1. - p[2] * p[2]) * (1. - p[0] * p[0]));
354 break;
355 default:
357 }
358 }
359 }
360
361 template <int dim>
363 PillowFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
364 {
365 Tensor<1, dim> result;
366 switch (dim)
367 {
368 case 1:
369 result[0] = -2. * p[0];
370 break;
371 case 2:
372 result[0] = -2. * p[0] * (1. - p[1] * p[1]);
373 result[1] = -2. * p[1] * (1. - p[0] * p[0]);
374 break;
375 case 3:
376 result[0] = -2. * p[0] * (1. - p[1] * p[1]) * (1. - p[2] * p[2]);
377 result[1] = -2. * p[1] * (1. - p[0] * p[0]) * (1. - p[2] * p[2]);
378 result[2] = -2. * p[2] * (1. - p[0] * p[0]) * (1. - p[1] * p[1]);
379 break;
380 default:
382 }
383 return result;
384 }
385
386 template <int dim>
387 void
389 std::vector<Tensor<1, dim>> &gradients,
390 const unsigned int) const
391 {
392 Assert(gradients.size() == points.size(),
393 ExcDimensionMismatch(gradients.size(), points.size()));
394
395 for (unsigned int i = 0; i < points.size(); ++i)
396 {
397 const Point<dim> &p = points[i];
398 switch (dim)
399 {
400 case 1:
401 gradients[i][0] = -2. * p[0];
402 break;
403 case 2:
404 gradients[i][0] = -2. * p[0] * (1. - p[1] * p[1]);
405 gradients[i][1] = -2. * p[1] * (1. - p[0] * p[0]);
406 break;
407 case 3:
408 gradients[i][0] =
409 -2. * p[0] * (1. - p[1] * p[1]) * (1. - p[2] * p[2]);
410 gradients[i][1] =
411 -2. * p[1] * (1. - p[0] * p[0]) * (1. - p[2] * p[2]);
412 gradients[i][2] =
413 -2. * p[2] * (1. - p[0] * p[0]) * (1. - p[1] * p[1]);
414 break;
415 default:
417 }
418 }
419 }
420
421 //--------------------------------------------------------------------
422
423 template <int dim>
424 CosineFunction<dim>::CosineFunction(const unsigned int n_components)
425 : Function<dim>(n_components)
426 {}
427
428
429
430 template <int dim>
431 double
432 CosineFunction<dim>::value(const Point<dim> &p, const unsigned int) const
433 {
434 switch (dim)
435 {
436 case 1:
437 return std::cos(numbers::PI_2 * p[0]);
438 case 2:
439 return std::cos(numbers::PI_2 * p[0]) *
440 std::cos(numbers::PI_2 * p[1]);
441 case 3:
442 return std::cos(numbers::PI_2 * p[0]) *
443 std::cos(numbers::PI_2 * p[1]) *
444 std::cos(numbers::PI_2 * p[2]);
445 default:
447 }
448 return 0.;
449 }
450
451 template <int dim>
452 void
453 CosineFunction<dim>::value_list(const std::vector<Point<dim>> &points,
454 std::vector<double> &values,
455 const unsigned int) const
456 {
457 Assert(values.size() == points.size(),
458 ExcDimensionMismatch(values.size(), points.size()));
459
460 for (unsigned int i = 0; i < points.size(); ++i)
461 values[i] = value(points[i]);
462 }
463
464
465 template <int dim>
466 void
468 const std::vector<Point<dim>> &points,
469 std::vector<Vector<double>> &values) const
470 {
471 Assert(values.size() == points.size(),
472 ExcDimensionMismatch(values.size(), points.size()));
473
474 for (unsigned int i = 0; i < points.size(); ++i)
475 {
476 const double v = value(points[i]);
477 for (unsigned int k = 0; k < values[i].size(); ++k)
478 values[i](k) = v;
479 }
480 }
481
482
483 template <int dim>
484 double
485 CosineFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
486 {
487 switch (dim)
488 {
489 case 1:
490 return -numbers::PI_2 * numbers::PI_2 *
491 std::cos(numbers::PI_2 * p[0]);
492 case 2:
493 return -2 * numbers::PI_2 * numbers::PI_2 *
494 std::cos(numbers::PI_2 * p[0]) *
495 std::cos(numbers::PI_2 * p[1]);
496 case 3:
497 return -3 * numbers::PI_2 * numbers::PI_2 *
498 std::cos(numbers::PI_2 * p[0]) *
499 std::cos(numbers::PI_2 * p[1]) *
500 std::cos(numbers::PI_2 * p[2]);
501 default:
503 }
504 return 0.;
505 }
506
507 template <int dim>
508 void
510 std::vector<double> &values,
511 const unsigned int) const
512 {
513 Assert(values.size() == points.size(),
514 ExcDimensionMismatch(values.size(), points.size()));
515
516 for (unsigned int i = 0; i < points.size(); ++i)
517 values[i] = laplacian(points[i]);
518 }
519
520 template <int dim>
522 CosineFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
523 {
524 Tensor<1, dim> result;
525 switch (dim)
526 {
527 case 1:
528 result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]);
529 break;
530 case 2:
531 result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
532 std::cos(numbers::PI_2 * p[1]);
533 result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
534 std::sin(numbers::PI_2 * p[1]);
535 break;
536 case 3:
537 result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
538 std::cos(numbers::PI_2 * p[1]) *
539 std::cos(numbers::PI_2 * p[2]);
540 result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
541 std::sin(numbers::PI_2 * p[1]) *
542 std::cos(numbers::PI_2 * p[2]);
543 result[2] = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
544 std::cos(numbers::PI_2 * p[1]) *
545 std::sin(numbers::PI_2 * p[2]);
546 break;
547 default:
549 }
550 return result;
551 }
552
553 template <int dim>
554 void
556 std::vector<Tensor<1, dim>> &gradients,
557 const unsigned int) const
558 {
559 Assert(gradients.size() == points.size(),
560 ExcDimensionMismatch(gradients.size(), points.size()));
561
562 for (unsigned int i = 0; i < points.size(); ++i)
563 {
564 const Point<dim> &p = points[i];
565 switch (dim)
566 {
567 case 1:
568 gradients[i][0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]);
569 break;
570 case 2:
571 gradients[i][0] = -numbers::PI_2 *
572 std::sin(numbers::PI_2 * p[0]) *
573 std::cos(numbers::PI_2 * p[1]);
574 gradients[i][1] = -numbers::PI_2 *
575 std::cos(numbers::PI_2 * p[0]) *
576 std::sin(numbers::PI_2 * p[1]);
577 break;
578 case 3:
579 gradients[i][0] =
582 gradients[i][1] =
585 gradients[i][2] =
588 break;
589 default:
591 }
592 }
593 }
594
595 template <int dim>
597 CosineFunction<dim>::hessian(const Point<dim> &p, const unsigned int) const
598 {
599 const double pi2 = numbers::PI_2 * numbers::PI_2;
600
602 switch (dim)
603 {
604 case 1:
605 result[0][0] = -pi2 * std::cos(numbers::PI_2 * p[0]);
606 break;
607 case 2:
608 {
609 const double coco = -pi2 * std::cos(numbers::PI_2 * p[0]) *
610 std::cos(numbers::PI_2 * p[1]);
611 const double sisi = pi2 * std::sin(numbers::PI_2 * p[0]) *
612 std::sin(numbers::PI_2 * p[1]);
613 result[0][0] = coco;
614 result[1][1] = coco;
615 // for SymmetricTensor we assign [ij] and [ji] simultaneously:
616 result[0][1] = sisi;
617 }
618 break;
619 case 3:
620 {
621 const double cococo = -pi2 * std::cos(numbers::PI_2 * p[0]) *
622 std::cos(numbers::PI_2 * p[1]) *
623 std::cos(numbers::PI_2 * p[2]);
624 const double sisico = pi2 * std::sin(numbers::PI_2 * p[0]) *
625 std::sin(numbers::PI_2 * p[1]) *
626 std::cos(numbers::PI_2 * p[2]);
627 const double sicosi = pi2 * std::sin(numbers::PI_2 * p[0]) *
628 std::cos(numbers::PI_2 * p[1]) *
629 std::sin(numbers::PI_2 * p[2]);
630 const double cosisi = pi2 * std::cos(numbers::PI_2 * p[0]) *
631 std::sin(numbers::PI_2 * p[1]) *
632 std::sin(numbers::PI_2 * p[2]);
633
634 result[0][0] = cococo;
635 result[1][1] = cococo;
636 result[2][2] = cococo;
637 // for SymmetricTensor we assign [ij] and [ji] simultaneously:
638 result[0][1] = sisico;
639 result[0][2] = sicosi;
640 result[1][2] = cosisi;
641 }
642 break;
643 default:
645 }
646 return result;
647 }
648
649 template <int dim>
650 void
652 const std::vector<Point<dim>> &points,
653 std::vector<SymmetricTensor<2, dim>> &hessians,
654 const unsigned int) const
655 {
656 Assert(hessians.size() == points.size(),
657 ExcDimensionMismatch(hessians.size(), points.size()));
658
659 const double pi2 = numbers::PI_2 * numbers::PI_2;
660
661 for (unsigned int i = 0; i < points.size(); ++i)
662 {
663 const Point<dim> &p = points[i];
664 switch (dim)
665 {
666 case 1:
667 hessians[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p[0]);
668 break;
669 case 2:
670 {
671 const double coco = -pi2 * std::cos(numbers::PI_2 * p[0]) *
672 std::cos(numbers::PI_2 * p[1]);
673 const double sisi = pi2 * std::sin(numbers::PI_2 * p[0]) *
674 std::sin(numbers::PI_2 * p[1]);
675 hessians[i][0][0] = coco;
676 hessians[i][1][1] = coco;
677 // for SymmetricTensor we assign [ij] and [ji] simultaneously:
678 hessians[i][0][1] = sisi;
679 }
680 break;
681 case 3:
682 {
683 const double cococo = -pi2 * std::cos(numbers::PI_2 * p[0]) *
684 std::cos(numbers::PI_2 * p[1]) *
685 std::cos(numbers::PI_2 * p[2]);
686 const double sisico = pi2 * std::sin(numbers::PI_2 * p[0]) *
687 std::sin(numbers::PI_2 * p[1]) *
688 std::cos(numbers::PI_2 * p[2]);
689 const double sicosi = pi2 * std::sin(numbers::PI_2 * p[0]) *
690 std::cos(numbers::PI_2 * p[1]) *
691 std::sin(numbers::PI_2 * p[2]);
692 const double cosisi = pi2 * std::cos(numbers::PI_2 * p[0]) *
693 std::sin(numbers::PI_2 * p[1]) *
694 std::sin(numbers::PI_2 * p[2]);
695
696 hessians[i][0][0] = cococo;
697 hessians[i][1][1] = cococo;
698 hessians[i][2][2] = cococo;
699 // for SymmetricTensor we assign [ij] and [ji] simultaneously:
700 hessians[i][0][1] = sisico;
701 hessians[i][0][2] = sicosi;
702 hessians[i][1][2] = cosisi;
703 }
704 break;
705 default:
707 }
708 }
709 }
710
711 //--------------------------------------------------------------------
712
713 template <int dim>
717
718
719 template <int dim>
720 double
722 const unsigned int d) const
723 {
724 AssertIndexRange(d, dim);
725 const unsigned int d1 = (d + 1) % dim;
726 const unsigned int d2 = (d + 2) % dim;
727 switch (dim)
728 {
729 case 1:
730 return (-numbers::PI_2 * std::sin(numbers::PI_2 * p[0]));
731 case 2:
732 return (-numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) *
733 std::cos(numbers::PI_2 * p[d1]));
734 case 3:
735 return (-numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) *
736 std::cos(numbers::PI_2 * p[d1]) *
737 std::cos(numbers::PI_2 * p[d2]));
738 default:
740 }
741 return 0.;
742 }
743
744
745 template <int dim>
746 void
748 Vector<double> &result) const
749 {
750 AssertDimension(result.size(), dim);
751 switch (dim)
752 {
753 case 1:
754 result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]);
755 break;
756 case 2:
757 result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
758 std::cos(numbers::PI_2 * p[1]);
759 result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
760 std::sin(numbers::PI_2 * p[1]);
761 break;
762 case 3:
763 result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
764 std::cos(numbers::PI_2 * p[1]) *
765 std::cos(numbers::PI_2 * p[2]);
766 result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
767 std::sin(numbers::PI_2 * p[1]) *
768 std::cos(numbers::PI_2 * p[2]);
769 result(2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
770 std::cos(numbers::PI_2 * p[1]) *
771 std::sin(numbers::PI_2 * p[2]);
772 break;
773 default:
775 }
776 }
777
778
779 template <int dim>
780 void
782 std::vector<double> &values,
783 const unsigned int d) const
784 {
785 Assert(values.size() == points.size(),
786 ExcDimensionMismatch(values.size(), points.size()));
787 AssertIndexRange(d, dim);
788 const unsigned int d1 = (d + 1) % dim;
789 const unsigned int d2 = (d + 2) % dim;
790
791 for (unsigned int i = 0; i < points.size(); ++i)
792 {
793 const Point<dim> &p = points[i];
794 switch (dim)
795 {
796 case 1:
797 values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[d]);
798 break;
799 case 2:
800 values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) *
801 std::cos(numbers::PI_2 * p[d1]);
802 break;
803 case 3:
804 values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) *
805 std::cos(numbers::PI_2 * p[d1]) *
806 std::cos(numbers::PI_2 * p[d2]);
807 break;
808 default:
810 }
811 }
812 }
813
814
815 template <int dim>
816 void
818 const std::vector<Point<dim>> &points,
819 std::vector<Vector<double>> &values) const
820 {
821 Assert(values.size() == points.size(),
822 ExcDimensionMismatch(values.size(), points.size()));
823
824 for (unsigned int i = 0; i < points.size(); ++i)
825 {
826 const Point<dim> &p = points[i];
827 switch (dim)
828 {
829 case 1:
830 values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]);
831 break;
832 case 2:
833 values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
834 std::cos(numbers::PI_2 * p[1]);
835 values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
836 std::sin(numbers::PI_2 * p[1]);
837 break;
838 case 3:
839 values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
840 std::cos(numbers::PI_2 * p[1]) *
841 std::cos(numbers::PI_2 * p[2]);
842 values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
843 std::sin(numbers::PI_2 * p[1]) *
844 std::cos(numbers::PI_2 * p[2]);
845 values[i](2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
846 std::cos(numbers::PI_2 * p[1]) *
847 std::sin(numbers::PI_2 * p[2]);
848 break;
849 default:
851 }
852 }
853 }
854
855
856 template <int dim>
857 double
859 const unsigned int d) const
860 {
861 return -numbers::PI_2 * numbers::PI_2 * value(p, d);
862 }
863
864
865 template <int dim>
868 const unsigned int d) const
869 {
870 AssertIndexRange(d, dim);
871 const unsigned int d1 = (d + 1) % dim;
872 const unsigned int d2 = (d + 2) % dim;
873 const double pi2 = numbers::PI_2 * numbers::PI_2;
874
875 Tensor<1, dim> result;
876 switch (dim)
877 {
878 case 1:
879 result[0] = -pi2 * std::cos(numbers::PI_2 * p[0]);
880 break;
881 case 2:
882 result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) *
883 std::cos(numbers::PI_2 * p[d1]);
884 result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) *
885 std::sin(numbers::PI_2 * p[d1]);
886 break;
887 case 3:
888 result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) *
889 std::cos(numbers::PI_2 * p[d1]) *
890 std::cos(numbers::PI_2 * p[d2]);
891 result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) *
892 std::sin(numbers::PI_2 * p[d1]) *
893 std::cos(numbers::PI_2 * p[d2]);
894 result[d2] = pi2 * std::sin(numbers::PI_2 * p[d]) *
895 std::cos(numbers::PI_2 * p[d1]) *
896 std::sin(numbers::PI_2 * p[d2]);
897 break;
898 default:
900 }
901 return result;
902 }
903
904
905 template <int dim>
906 void
908 std::vector<Tensor<1, dim>> &gradients,
909 const unsigned int d) const
910 {
911 AssertIndexRange(d, dim);
912 const unsigned int d1 = (d + 1) % dim;
913 const unsigned int d2 = (d + 2) % dim;
914 const double pi2 = numbers::PI_2 * numbers::PI_2;
915
916 Assert(gradients.size() == points.size(),
917 ExcDimensionMismatch(gradients.size(), points.size()));
918 for (unsigned int i = 0; i < points.size(); ++i)
919 {
920 const Point<dim> &p = points[i];
921 Tensor<1, dim> &result = gradients[i];
922
923 switch (dim)
924 {
925 case 1:
926 result[0] = -pi2 * std::cos(numbers::PI_2 * p[0]);
927 break;
928 case 2:
929 result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) *
930 std::cos(numbers::PI_2 * p[d1]);
931 result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) *
932 std::sin(numbers::PI_2 * p[d1]);
933 break;
934 case 3:
935 result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) *
936 std::cos(numbers::PI_2 * p[d1]) *
937 std::cos(numbers::PI_2 * p[d2]);
938 result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) *
939 std::sin(numbers::PI_2 * p[d1]) *
940 std::cos(numbers::PI_2 * p[d2]);
941 result[d2] = pi2 * std::sin(numbers::PI_2 * p[d]) *
942 std::cos(numbers::PI_2 * p[d1]) *
943 std::sin(numbers::PI_2 * p[d2]);
944 break;
945 default:
947 }
948 }
949 }
950
951
952 template <int dim>
953 void
955 const std::vector<Point<dim>> &points,
956 std::vector<std::vector<Tensor<1, dim>>> &gradients) const
957 {
958 AssertVectorVectorDimension(gradients, points.size(), dim);
959 const double pi2 = numbers::PI_2 * numbers::PI_2;
960
961 for (unsigned int i = 0; i < points.size(); ++i)
962 {
963 const Point<dim> &p = points[i];
964 switch (dim)
965 {
966 case 1:
967 gradients[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p[0]);
968 break;
969 case 2:
970 {
971 const double coco = -pi2 * std::cos(numbers::PI_2 * p[0]) *
972 std::cos(numbers::PI_2 * p[1]);
973 const double sisi = pi2 * std::sin(numbers::PI_2 * p[0]) *
974 std::sin(numbers::PI_2 * p[1]);
975 gradients[i][0][0] = coco;
976 gradients[i][1][1] = coco;
977 gradients[i][0][1] = sisi;
978 gradients[i][1][0] = sisi;
979 }
980 break;
981 case 3:
982 {
983 const double cococo = -pi2 * std::cos(numbers::PI_2 * p[0]) *
984 std::cos(numbers::PI_2 * p[1]) *
985 std::cos(numbers::PI_2 * p[2]);
986 const double sisico = pi2 * std::sin(numbers::PI_2 * p[0]) *
987 std::sin(numbers::PI_2 * p[1]) *
988 std::cos(numbers::PI_2 * p[2]);
989 const double sicosi = pi2 * std::sin(numbers::PI_2 * p[0]) *
990 std::cos(numbers::PI_2 * p[1]) *
991 std::sin(numbers::PI_2 * p[2]);
992 const double cosisi = pi2 * std::cos(numbers::PI_2 * p[0]) *
993 std::sin(numbers::PI_2 * p[1]) *
994 std::sin(numbers::PI_2 * p[2]);
995
996 gradients[i][0][0] = cococo;
997 gradients[i][1][1] = cococo;
998 gradients[i][2][2] = cococo;
999 gradients[i][0][1] = sisico;
1000 gradients[i][1][0] = sisico;
1001 gradients[i][0][2] = sicosi;
1002 gradients[i][2][0] = sicosi;
1003 gradients[i][1][2] = cosisi;
1004 gradients[i][2][1] = cosisi;
1005 }
1006 break;
1007 default:
1009 }
1010 }
1011 }
1012
1013
1014 //--------------------------------------------------------------------
1015
1016 template <int dim>
1017 double
1018 ExpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1019 {
1020 switch (dim)
1021 {
1022 case 1:
1023 return std::exp(p[0]);
1024 case 2:
1025 return std::exp(p[0]) * std::exp(p[1]);
1026 case 3:
1027 return std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]);
1028 default:
1030 }
1031 return 0.;
1032 }
1033
1034 template <int dim>
1035 void
1036 ExpFunction<dim>::value_list(const std::vector<Point<dim>> &points,
1037 std::vector<double> &values,
1038 const unsigned int) const
1039 {
1040 Assert(values.size() == points.size(),
1041 ExcDimensionMismatch(values.size(), points.size()));
1042
1043 for (unsigned int i = 0; i < points.size(); ++i)
1044 {
1045 const Point<dim> &p = points[i];
1046 switch (dim)
1047 {
1048 case 1:
1049 values[i] = std::exp(p[0]);
1050 break;
1051 case 2:
1052 values[i] = std::exp(p[0]) * std::exp(p[1]);
1053 break;
1054 case 3:
1055 values[i] = std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]);
1056 break;
1057 default:
1059 }
1060 }
1061 }
1062
1063 template <int dim>
1064 double
1065 ExpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1066 {
1067 switch (dim)
1068 {
1069 case 1:
1070 return std::exp(p[0]);
1071 case 2:
1072 return 2 * std::exp(p[0]) * std::exp(p[1]);
1073 case 3:
1074 return 3 * std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]);
1075 default:
1077 }
1078 return 0.;
1079 }
1080
1081 template <int dim>
1082 void
1084 std::vector<double> &values,
1085 const unsigned int) const
1086 {
1087 Assert(values.size() == points.size(),
1088 ExcDimensionMismatch(values.size(), points.size()));
1089
1090 for (unsigned int i = 0; i < points.size(); ++i)
1091 {
1092 const Point<dim> &p = points[i];
1093 switch (dim)
1094 {
1095 case 1:
1096 values[i] = std::exp(p[0]);
1097 break;
1098 case 2:
1099 values[i] = 2 * std::exp(p[0]) * std::exp(p[1]);
1100 break;
1101 case 3:
1102 values[i] = 3 * std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]);
1103 break;
1104 default:
1106 }
1107 }
1108 }
1109
1110 template <int dim>
1112 ExpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1113 {
1114 Tensor<1, dim> result;
1115 switch (dim)
1116 {
1117 case 1:
1118 result[0] = std::exp(p[0]);
1119 break;
1120 case 2:
1121 result[0] = std::exp(p[0]) * std::exp(p[1]);
1122 result[1] = result[0];
1123 break;
1124 case 3:
1125 result[0] = std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]);
1126 result[1] = result[0];
1127 result[2] = result[0];
1128 break;
1129 default:
1131 }
1132 return result;
1133 }
1134
1135 template <int dim>
1136 void
1138 std::vector<Tensor<1, dim>> &gradients,
1139 const unsigned int) const
1140 {
1141 Assert(gradients.size() == points.size(),
1142 ExcDimensionMismatch(gradients.size(), points.size()));
1143
1144 for (unsigned int i = 0; i < points.size(); ++i)
1145 {
1146 const Point<dim> &p = points[i];
1147 switch (dim)
1148 {
1149 case 1:
1150 gradients[i][0] = std::exp(p[0]);
1151 break;
1152 case 2:
1153 gradients[i][0] = std::exp(p[0]) * std::exp(p[1]);
1154 gradients[i][1] = gradients[i][0];
1155 break;
1156 case 3:
1157 gradients[i][0] =
1158 std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]);
1159 gradients[i][1] = gradients[i][0];
1160 gradients[i][2] = gradients[i][0];
1161 break;
1162 default:
1164 }
1165 }
1166 }
1167
1168 //--------------------------------------------------------------------
1169
1170
1171 double
1172 LSingularityFunction::value(const Point<2> &p, const unsigned int) const
1173 {
1174 const double x = p[0];
1175 const double y = p[1];
1176
1177 if ((x >= 0) && (y >= 0))
1178 return 0.;
1179
1180 const double phi = std::atan2(y, -x) + numbers::PI;
1181 const double r_squared = x * x + y * y;
1182
1183 return std::cbrt(r_squared) * std::sin(2. / 3. * phi);
1184 }
1185
1186
1187
1188 void
1189 LSingularityFunction::value_list(const std::vector<Point<2>> &points,
1190 std::vector<double> &values,
1191 const unsigned int) const
1192 {
1193 Assert(values.size() == points.size(),
1194 ExcDimensionMismatch(values.size(), points.size()));
1195
1196 for (unsigned int i = 0; i < points.size(); ++i)
1197 {
1198 const double x = points[i][0];
1199 const double y = points[i][1];
1200
1201 if ((x >= 0) && (y >= 0))
1202 values[i] = 0.;
1203 else
1204 {
1205 const double phi = std::atan2(y, -x) + numbers::PI;
1206 const double r_squared = x * x + y * y;
1207
1208 values[i] = std::cbrt(r_squared) * std::sin(2. / 3. * phi);
1209 }
1210 }
1211 }
1212
1213
1214
1215 void
1217 const std::vector<Point<2>> &points,
1218 std::vector<Vector<double>> &values) const
1219 {
1220 Assert(values.size() == points.size(),
1221 ExcDimensionMismatch(values.size(), points.size()));
1222
1223 for (unsigned int i = 0; i < points.size(); ++i)
1224 {
1225 Assert(values[i].size() == 1,
1226 ExcDimensionMismatch(values[i].size(), 1));
1227 const double x = points[i][0];
1228 const double y = points[i][1];
1229
1230 if ((x >= 0) && (y >= 0))
1231 values[i](0) = 0.;
1232 else
1233 {
1234 const double phi = std::atan2(y, -x) + numbers::PI;
1235 const double r_squared = x * x + y * y;
1236
1237 values[i](0) = std::cbrt(r_squared) * std::sin(2. / 3. * phi);
1238 }
1239 }
1240 }
1241
1242
1243
1244 double
1245 LSingularityFunction::laplacian(const Point<2> &, const unsigned int) const
1246 {
1247 // Not a bug but exactly how the function is defined:
1248 return 0.;
1249 }
1250
1251
1252
1253 void
1255 std::vector<double> &values,
1256 const unsigned int) const
1257 {
1258 Assert(values.size() == points.size(),
1259 ExcDimensionMismatch(values.size(), points.size()));
1260
1261 for (unsigned int i = 0; i < points.size(); ++i)
1262 values[i] = 0.;
1263 }
1264
1265
1266
1268 LSingularityFunction::gradient(const Point<2> &p, const unsigned int) const
1269 {
1270 const double x = p[0];
1271 const double y = p[1];
1272 const double phi = std::atan2(y, -x) + numbers::PI;
1273 const double r43 = std::pow(x * x + y * y, 2. / 3.);
1274
1275 Tensor<1, 2> result;
1276 result[0] = 2. / 3. *
1277 (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) /
1278 r43;
1279 result[1] = 2. / 3. *
1280 (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) /
1281 r43;
1282 return result;
1283 }
1284
1285
1286
1287 void
1289 std::vector<Tensor<1, 2>> &gradients,
1290 const unsigned int) const
1291 {
1292 Assert(gradients.size() == points.size(),
1293 ExcDimensionMismatch(gradients.size(), points.size()));
1294
1295 for (unsigned int i = 0; i < points.size(); ++i)
1296 {
1297 const Point<2> &p = points[i];
1298 const double x = p[0];
1299 const double y = p[1];
1300 const double phi = std::atan2(y, -x) + numbers::PI;
1301 const double r43 = std::pow(x * x + y * y, 2. / 3.);
1302
1303 gradients[i][0] =
1304 2. / 3. *
1305 (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1306 gradients[i][1] =
1307 2. / 3. *
1308 (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1309 }
1310 }
1311
1312
1313
1314 void
1316 const std::vector<Point<2>> &points,
1317 std::vector<std::vector<Tensor<1, 2>>> &gradients) const
1318 {
1319 Assert(gradients.size() == points.size(),
1320 ExcDimensionMismatch(gradients.size(), points.size()));
1321
1322 for (unsigned int i = 0; i < points.size(); ++i)
1323 {
1324 Assert(gradients[i].size() == 1,
1325 ExcDimensionMismatch(gradients[i].size(), 1));
1326 const Point<2> &p = points[i];
1327 const double x = p[0];
1328 const double y = p[1];
1329 const double phi = std::atan2(y, -x) + numbers::PI;
1330 const double r43 = std::pow(x * x + y * y, 2. / 3.);
1331
1332 gradients[i][0][0] =
1333 2. / 3. *
1334 (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1335 gradients[i][0][1] =
1336 2. / 3. *
1337 (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1338 }
1339 }
1340
1341 //--------------------------------------------------------------------
1342
1346
1347
1348
1349 double
1350 LSingularityGradFunction::value(const Point<2> &p, const unsigned int d) const
1351 {
1352 AssertIndexRange(d, 2);
1353
1354 const double x = p[0];
1355 const double y = p[1];
1356 const double phi = std::atan2(y, -x) + numbers::PI;
1357 const double r43 = std::pow(x * x + y * y, 2. / 3.);
1358
1359 return 2. / 3. *
1360 (std::sin(2. / 3. * phi) * p[d] +
1361 (d == 0 ? (std::cos(2. / 3. * phi) * y) :
1362 (-std::cos(2. / 3. * phi) * x))) /
1363 r43;
1364 }
1365
1366
1367 void
1369 std::vector<double> &values,
1370 const unsigned int d) const
1371 {
1372 AssertIndexRange(d, 2);
1373 AssertDimension(values.size(), points.size());
1374
1375 for (unsigned int i = 0; i < points.size(); ++i)
1376 {
1377 const Point<2> &p = points[i];
1378 const double x = p[0];
1379 const double y = p[1];
1380 const double phi = std::atan2(y, -x) + numbers::PI;
1381 const double r43 = std::pow(x * x + y * y, 2. / 3.);
1382
1383 values[i] = 2. / 3. *
1384 (std::sin(2. / 3. * phi) * p[d] +
1385 (d == 0 ? (std::cos(2. / 3. * phi) * y) :
1386 (-std::cos(2. / 3. * phi) * x))) /
1387 r43;
1388 }
1389 }
1390
1391
1392 void
1394 const std::vector<Point<2>> &points,
1395 std::vector<Vector<double>> &values) const
1396 {
1397 Assert(values.size() == points.size(),
1398 ExcDimensionMismatch(values.size(), points.size()));
1399
1400 for (unsigned int i = 0; i < points.size(); ++i)
1401 {
1402 AssertDimension(values[i].size(), 2);
1403 const Point<2> &p = points[i];
1404 const double x = p[0];
1405 const double y = p[1];
1406 const double phi = std::atan2(y, -x) + numbers::PI;
1407 const double r43 = std::pow(x * x + y * y, 2. / 3.);
1408
1409 values[i](0) =
1410 2. / 3. *
1411 (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1412 values[i](1) =
1413 2. / 3. *
1414 (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1415 }
1416 }
1417
1418
1419 double
1421 const unsigned int) const
1422 {
1423 return 0.;
1424 }
1425
1426
1427 void
1429 std::vector<double> &values,
1430 const unsigned int) const
1431 {
1432 Assert(values.size() == points.size(),
1433 ExcDimensionMismatch(values.size(), points.size()));
1434
1435 for (unsigned int i = 0; i < points.size(); ++i)
1436 values[i] = 0.;
1437 }
1438
1439
1440
1443 const unsigned int /*component*/) const
1444 {
1446 return {};
1447 }
1448
1449
1450 void
1452 const std::vector<Point<2>> & /*points*/,
1453 std::vector<Tensor<1, 2>> & /*gradients*/,
1454 const unsigned int /*component*/) const
1455 {
1457 }
1458
1459
1460 void
1462 const std::vector<Point<2>> & /*points*/,
1463 std::vector<std::vector<Tensor<1, 2>>> & /*gradients*/) const
1464 {
1466 }
1467
1468 //--------------------------------------------------------------------
1469
1470 template <int dim>
1471 double
1473 const unsigned int) const
1474 {
1475 const double x = p[0];
1476 const double y = p[1];
1477
1478 const double phi = std::atan2(x, y) + numbers::PI;
1479 const double r_squared = x * x + y * y;
1480
1481 return std::pow(r_squared, .25) * std::sin(.5 * phi);
1482 }
1483
1484
1485 template <int dim>
1486 void
1488 const std::vector<Point<dim>> &points,
1489 std::vector<double> &values,
1490 const unsigned int) const
1491 {
1492 Assert(values.size() == points.size(),
1493 ExcDimensionMismatch(values.size(), points.size()));
1494
1495 for (unsigned int i = 0; i < points.size(); ++i)
1496 {
1497 const double x = points[i][0];
1498 const double y = points[i][1];
1499
1500 const double phi = std::atan2(x, y) + numbers::PI;
1501 const double r_squared = x * x + y * y;
1502
1503 values[i] = std::pow(r_squared, .25) * std::sin(.5 * phi);
1504 }
1505 }
1506
1507
1508 template <int dim>
1509 void
1511 const std::vector<Point<dim>> &points,
1512 std::vector<Vector<double>> &values) const
1513 {
1514 Assert(values.size() == points.size(),
1515 ExcDimensionMismatch(values.size(), points.size()));
1516
1517 for (unsigned int i = 0; i < points.size(); ++i)
1518 {
1519 Assert(values[i].size() == 1,
1520 ExcDimensionMismatch(values[i].size(), 1));
1521
1522 const double x = points[i][0];
1523 const double y = points[i][1];
1524
1525 const double phi = std::atan2(x, y) + numbers::PI;
1526 const double r_squared = x * x + y * y;
1527
1528 values[i](0) = std::pow(r_squared, .25) * std::sin(.5 * phi);
1529 }
1530 }
1531
1532
1533 template <int dim>
1534 double
1536 const unsigned int) const
1537 {
1538 return 0.;
1539 }
1540
1541
1542 template <int dim>
1543 void
1545 const std::vector<Point<dim>> &points,
1546 std::vector<double> &values,
1547 const unsigned int) const
1548 {
1549 Assert(values.size() == points.size(),
1550 ExcDimensionMismatch(values.size(), points.size()));
1551
1552 for (unsigned int i = 0; i < points.size(); ++i)
1553 values[i] = 0.;
1554 }
1555
1556
1557 template <int dim>
1560 const unsigned int) const
1561 {
1562 const double x = p[0];
1563 const double y = p[1];
1564 const double phi = std::atan2(x, y) + numbers::PI;
1565 const double r64 = std::pow(x * x + y * y, 3. / 4.);
1566
1567 Tensor<1, dim> result;
1568 result[0] = 1. / 2. *
1569 (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) /
1570 r64;
1571 result[1] = 1. / 2. *
1572 (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) /
1573 r64;
1574 return result;
1575 }
1576
1577
1578 template <int dim>
1579 void
1581 const std::vector<Point<dim>> &points,
1582 std::vector<Tensor<1, dim>> &gradients,
1583 const unsigned int) const
1584 {
1585 Assert(gradients.size() == points.size(),
1586 ExcDimensionMismatch(gradients.size(), points.size()));
1587
1588 for (unsigned int i = 0; i < points.size(); ++i)
1589 {
1590 const Point<dim> &p = points[i];
1591 const double x = p[0];
1592 const double y = p[1];
1593 const double phi = std::atan2(x, y) + numbers::PI;
1594 const double r64 = std::pow(x * x + y * y, 3. / 4.);
1595
1596 gradients[i][0] =
1597 1. / 2. *
1598 (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) / r64;
1599 gradients[i][1] =
1600 1. / 2. *
1601 (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) / r64;
1602 for (unsigned int d = 2; d < dim; ++d)
1603 gradients[i][d] = 0.;
1604 }
1605 }
1606
1607 template <int dim>
1608 void
1610 const std::vector<Point<dim>> &points,
1611 std::vector<std::vector<Tensor<1, dim>>> &gradients) const
1612 {
1613 Assert(gradients.size() == points.size(),
1614 ExcDimensionMismatch(gradients.size(), points.size()));
1615
1616 for (unsigned int i = 0; i < points.size(); ++i)
1617 {
1618 Assert(gradients[i].size() == 1,
1619 ExcDimensionMismatch(gradients[i].size(), 1));
1620
1621 const Point<dim> &p = points[i];
1622 const double x = p[0];
1623 const double y = p[1];
1624 const double phi = std::atan2(x, y) + numbers::PI;
1625 const double r64 = std::pow(x * x + y * y, 3. / 4.);
1626
1627 gradients[i][0][0] =
1628 1. / 2. *
1629 (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) / r64;
1630 gradients[i][0][1] =
1631 1. / 2. *
1632 (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) / r64;
1633 for (unsigned int d = 2; d < dim; ++d)
1634 gradients[i][0][d] = 0.;
1635 }
1636 }
1637
1638 //--------------------------------------------------------------------
1639
1640
1641 double
1643 const unsigned int) const
1644 {
1645 const double x = p[0];
1646 const double y = p[1];
1647
1648 const double phi = std::atan2(x, y) + numbers::PI;
1649 const double r_squared = x * x + y * y;
1650
1651 return std::pow(r_squared, .125) * std::sin(.25 * phi);
1652 }
1653
1654
1655 void
1657 std::vector<double> &values,
1658 const unsigned int) const
1659 {
1660 Assert(values.size() == points.size(),
1661 ExcDimensionMismatch(values.size(), points.size()));
1662
1663 for (unsigned int i = 0; i < points.size(); ++i)
1664 {
1665 const double x = points[i][0];
1666 const double y = points[i][1];
1667
1668 const double phi = std::atan2(x, y) + numbers::PI;
1669 const double r_squared = x * x + y * y;
1670
1671 values[i] = std::pow(r_squared, .125) * std::sin(.25 * phi);
1672 }
1673 }
1674
1675
1676 void
1678 const std::vector<Point<2>> &points,
1679 std::vector<Vector<double>> &values) const
1680 {
1681 Assert(values.size() == points.size(),
1682 ExcDimensionMismatch(values.size(), points.size()));
1683
1684 for (unsigned int i = 0; i < points.size(); ++i)
1685 {
1686 Assert(values[i].size() == 1,
1687 ExcDimensionMismatch(values[i].size(), 1));
1688
1689 const double x = points[i][0];
1690 const double y = points[i][1];
1691
1692 const double phi = std::atan2(x, y) + numbers::PI;
1693 const double r_squared = x * x + y * y;
1694
1695 values[i](0) = std::pow(r_squared, .125) * std::sin(.25 * phi);
1696 }
1697 }
1698
1699
1700 double
1702 const unsigned int) const
1703 {
1704 return 0.;
1705 }
1706
1707
1708 void
1710 const std::vector<Point<2>> &points,
1711 std::vector<double> &values,
1712 const unsigned int) const
1713 {
1714 Assert(values.size() == points.size(),
1715 ExcDimensionMismatch(values.size(), points.size()));
1716
1717 for (unsigned int i = 0; i < points.size(); ++i)
1718 values[i] = 0.;
1719 }
1720
1721
1724 const unsigned int) const
1725 {
1726 const double x = p[0];
1727 const double y = p[1];
1728 const double phi = std::atan2(x, y) + numbers::PI;
1729 const double r78 = std::pow(x * x + y * y, 7. / 8.);
1730
1731
1732 Tensor<1, 2> result;
1733 result[0] = 1. / 4. *
1734 (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) /
1735 r78;
1736 result[1] = 1. / 4. *
1737 (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) /
1738 r78;
1739 return result;
1740 }
1741
1742
1743 void
1745 const std::vector<Point<2>> &points,
1746 std::vector<Tensor<1, 2>> &gradients,
1747 const unsigned int) const
1748 {
1749 Assert(gradients.size() == points.size(),
1750 ExcDimensionMismatch(gradients.size(), points.size()));
1751
1752 for (unsigned int i = 0; i < points.size(); ++i)
1753 {
1754 const Point<2> &p = points[i];
1755 const double x = p[0];
1756 const double y = p[1];
1757 const double phi = std::atan2(x, y) + numbers::PI;
1758 const double r78 = std::pow(x * x + y * y, 7. / 8.);
1759
1760 gradients[i][0] =
1761 1. / 4. *
1762 (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) / r78;
1763 gradients[i][1] =
1764 1. / 4. *
1765 (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) / r78;
1766 }
1767 }
1768
1769
1770 void
1772 const std::vector<Point<2>> &points,
1773 std::vector<std::vector<Tensor<1, 2>>> &gradients) const
1774 {
1775 Assert(gradients.size() == points.size(),
1776 ExcDimensionMismatch(gradients.size(), points.size()));
1777
1778 for (unsigned int i = 0; i < points.size(); ++i)
1779 {
1780 Assert(gradients[i].size() == 1,
1781 ExcDimensionMismatch(gradients[i].size(), 1));
1782
1783 const Point<2> &p = points[i];
1784 const double x = p[0];
1785 const double y = p[1];
1786 const double phi = std::atan2(x, y) + numbers::PI;
1787 const double r78 = std::pow(x * x + y * y, 7. / 8.);
1788
1789 gradients[i][0][0] =
1790 1. / 4. *
1791 (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) / r78;
1792 gradients[i][0][1] =
1793 1. / 4. *
1794 (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) / r78;
1795 }
1796 }
1797
1798 //--------------------------------------------------------------------
1799
1800 template <int dim>
1802 const double steepness)
1803 : direction(direction)
1804 , steepness(steepness)
1805 {
1806 switch (dim)
1807 {
1808 case 1:
1809 angle = 0;
1810 break;
1811 case 2:
1812 angle = std::atan2(direction[0], direction[1]);
1813 break;
1814 case 3:
1816 }
1817 sine = std::sin(angle);
1819 }
1820
1821
1822
1823 template <int dim>
1824 double
1825 JumpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1826 {
1827 const double x = steepness * (-cosine * p[0] + sine * p[1]);
1828 return -std::atan(x);
1829 }
1830
1831
1832
1833 template <int dim>
1834 void
1836 std::vector<double> &values,
1837 const unsigned int) const
1838 {
1839 Assert(values.size() == p.size(),
1840 ExcDimensionMismatch(values.size(), p.size()));
1841
1842 for (unsigned int i = 0; i < p.size(); ++i)
1843 {
1844 const double x = steepness * (-cosine * p[i][0] + sine * p[i][1]);
1845 values[i] = -std::atan(x);
1846 }
1847 }
1848
1849
1850 template <int dim>
1851 double
1852 JumpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1853 {
1854 const double x = steepness * (-cosine * p[0] + sine * p[1]);
1855 const double r = 1 + x * x;
1856 return 2 * steepness * steepness * x / (r * r);
1857 }
1858
1859
1860 template <int dim>
1861 void
1863 std::vector<double> &values,
1864 const unsigned int) const
1865 {
1866 Assert(values.size() == p.size(),
1867 ExcDimensionMismatch(values.size(), p.size()));
1868
1869 double f = 2 * steepness * steepness;
1870
1871 for (unsigned int i = 0; i < p.size(); ++i)
1872 {
1873 const double x = steepness * (-cosine * p[i][0] + sine * p[i][1]);
1874 const double r = 1 + x * x;
1875 values[i] = f * x / (r * r);
1876 }
1877 }
1878
1879
1880
1881 template <int dim>
1883 JumpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1884 {
1885 const double x = steepness * (-cosine * p[0] + sine * p[1]);
1886 const double r = -steepness * (1 + x * x);
1887 Tensor<1, dim> erg;
1888 erg[0] = cosine * r;
1889 erg[1] = sine * r;
1890 return erg;
1891 }
1892
1893
1894
1895 template <int dim>
1896 void
1898 std::vector<Tensor<1, dim>> &gradients,
1899 const unsigned int) const
1900 {
1901 Assert(gradients.size() == p.size(),
1902 ExcDimensionMismatch(gradients.size(), p.size()));
1903
1904 for (unsigned int i = 0; i < p.size(); ++i)
1905 {
1906 const double x = steepness * (cosine * p[i][0] + sine * p[i][1]);
1907 const double r = -steepness * (1 + x * x);
1908 gradients[i][0] = cosine * r;
1909 gradients[i][1] = sine * r;
1910 }
1911 }
1912
1913
1914
1915 template <int dim>
1916 std::size_t
1918 {
1919 // only simple data elements, so
1920 // use sizeof operator
1921 return sizeof(*this);
1922 }
1923
1924
1925
1926 /* ---------------------- FourierCosineFunction ----------------------- */
1927
1928
1929 template <int dim>
1931 const Tensor<1, dim> &fourier_coefficients)
1932 : Function<dim>(1)
1933 , fourier_coefficients(fourier_coefficients)
1934 {}
1935
1936
1937
1938 template <int dim>
1939 double
1941 const unsigned int component) const
1942 {
1943 (void)component;
1944 AssertIndexRange(component, 1);
1945 return std::cos(fourier_coefficients * p);
1946 }
1947
1948
1949
1950 template <int dim>
1953 const unsigned int component) const
1954 {
1955 (void)component;
1956 AssertIndexRange(component, 1);
1957 return -fourier_coefficients * std::sin(fourier_coefficients * p);
1958 }
1959
1960
1961
1962 template <int dim>
1963 double
1965 const unsigned int component) const
1966 {
1967 (void)component;
1968 AssertIndexRange(component, 1);
1969 return (fourier_coefficients * fourier_coefficients) *
1970 (-std::cos(fourier_coefficients * p));
1971 }
1972
1973
1974
1975 /* ---------------------- FourierSineFunction ----------------------- */
1976
1977
1978
1979 template <int dim>
1981 const Tensor<1, dim> &fourier_coefficients)
1982 : Function<dim>(1)
1983 , fourier_coefficients(fourier_coefficients)
1984 {}
1985
1986
1987
1988 template <int dim>
1989 double
1991 const unsigned int component) const
1992 {
1993 (void)component;
1994 AssertIndexRange(component, 1);
1995 return std::sin(fourier_coefficients * p);
1996 }
1997
1998
1999
2000 template <int dim>
2003 const unsigned int component) const
2004 {
2005 (void)component;
2006 AssertIndexRange(component, 1);
2007 return fourier_coefficients * std::cos(fourier_coefficients * p);
2008 }
2009
2010
2011
2012 template <int dim>
2013 double
2015 const unsigned int component) const
2016 {
2017 (void)component;
2018 AssertIndexRange(component, 1);
2019 return (fourier_coefficients * fourier_coefficients) *
2020 (-std::sin(fourier_coefficients * p));
2021 }
2022
2023
2024
2025 /* ---------------------- FourierSineSum ----------------------- */
2026
2027
2028
2029 template <int dim>
2031 const std::vector<Point<dim>> &fourier_coefficients,
2032 const std::vector<double> &weights)
2033 : Function<dim>(1)
2034 , fourier_coefficients(fourier_coefficients)
2035 , weights(weights)
2036 {
2037 Assert(fourier_coefficients.size() > 0, ExcZero());
2038 Assert(fourier_coefficients.size() == weights.size(),
2040 }
2041
2042
2043
2044 template <int dim>
2045 double
2047 const unsigned int component) const
2048 {
2049 (void)component;
2050 AssertIndexRange(component, 1);
2051
2052 const unsigned int n = weights.size();
2053 double sum = 0;
2054 for (unsigned int s = 0; s < n; ++s)
2055 sum += weights[s] * std::sin(fourier_coefficients[s] * p);
2056
2057 return sum;
2058 }
2059
2060
2061
2062 template <int dim>
2065 const unsigned int component) const
2066 {
2067 (void)component;
2068 AssertIndexRange(component, 1);
2069
2070 const unsigned int n = weights.size();
2071 Tensor<1, dim> sum;
2072 for (unsigned int s = 0; s < n; ++s)
2073 sum += fourier_coefficients[s] * std::cos(fourier_coefficients[s] * p);
2074
2075 return sum;
2076 }
2077
2078
2079
2080 template <int dim>
2081 double
2083 const unsigned int component) const
2084 {
2085 (void)component;
2086 AssertIndexRange(component, 1);
2087
2088 const unsigned int n = weights.size();
2089 double sum = 0;
2090 for (unsigned int s = 0; s < n; ++s)
2091 sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2092 std::sin(fourier_coefficients[s] * p);
2093
2094 return sum;
2095 }
2096
2097
2098
2099 /* ---------------------- FourierCosineSum ----------------------- */
2100
2101
2102
2103 template <int dim>
2105 const std::vector<Point<dim>> &fourier_coefficients,
2106 const std::vector<double> &weights)
2107 : Function<dim>(1)
2108 , fourier_coefficients(fourier_coefficients)
2109 , weights(weights)
2110 {
2111 Assert(fourier_coefficients.size() > 0, ExcZero());
2112 Assert(fourier_coefficients.size() == weights.size(),
2114 }
2115
2116
2117
2118 template <int dim>
2119 double
2121 const unsigned int component) const
2122 {
2123 (void)component;
2124 AssertIndexRange(component, 1);
2125
2126 const unsigned int n = weights.size();
2127 double sum = 0;
2128 for (unsigned int s = 0; s < n; ++s)
2129 sum += weights[s] * std::cos(fourier_coefficients[s] * p);
2130
2131 return sum;
2132 }
2133
2134
2135
2136 template <int dim>
2139 const unsigned int component) const
2140 {
2141 (void)component;
2142 AssertIndexRange(component, 1);
2143
2144 const unsigned int n = weights.size();
2145 Tensor<1, dim> sum;
2146 for (unsigned int s = 0; s < n; ++s)
2147 sum -= fourier_coefficients[s] * std::sin(fourier_coefficients[s] * p);
2148
2149 return sum;
2150 }
2151
2152
2153
2154 template <int dim>
2155 double
2157 const unsigned int component) const
2158 {
2159 (void)component;
2160 AssertIndexRange(component, 1);
2161
2162 const unsigned int n = weights.size();
2163 double sum = 0;
2164 for (unsigned int s = 0; s < n; ++s)
2165 sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2166 std::cos(fourier_coefficients[s] * p);
2167
2168 return sum;
2169 }
2170
2171
2172
2173 /* ---------------------- Monomial ----------------------- */
2174
2175
2176
2177 template <int dim, typename Number>
2179 const unsigned int n_components)
2180 : Function<dim, Number>(n_components)
2181 , exponents(exponents)
2182 {}
2183
2184
2185
2186 template <int dim, typename Number>
2187 Number
2189 const unsigned int component) const
2190 {
2191 (void)component;
2192 AssertIndexRange(component, this->n_components);
2193
2194 Number prod = 1;
2195 for (unsigned int s = 0; s < dim; ++s)
2196 {
2197 if (p[s] < 0)
2198 Assert(std::floor(exponents[s]) == exponents[s],
2199 ExcMessage("Exponentiation of a negative base number with "
2200 "a real exponent can't be performed."));
2201 prod *= std::pow(p[s], exponents[s]);
2202 }
2203 return prod;
2204 }
2205
2206
2207
2208 template <int dim, typename Number>
2209 void
2211 Vector<Number> &values) const
2212 {
2213 Assert(values.size() == this->n_components,
2214 ExcDimensionMismatch(values.size(), this->n_components));
2215
2216 for (unsigned int i = 0; i < values.size(); ++i)
2217 values(i) = Monomial<dim, Number>::value(p, i);
2218 }
2219
2220
2221
2222 template <int dim, typename Number>
2225 const unsigned int component) const
2226 {
2227 (void)component;
2228 AssertIndexRange(component, 1);
2229
2231 for (unsigned int d = 0; d < dim; ++d)
2232 {
2233 double prod = 1;
2234 for (unsigned int s = 0; s < dim; ++s)
2235 {
2236 if ((s == d) && (exponents[s] == 0) && (p[s] == 0))
2237 {
2238 prod = 0;
2239 break;
2240 }
2241 else
2242 {
2243 if (p[s] < 0)
2244 Assert(std::floor(exponents[s]) == exponents[s],
2245 ExcMessage(
2246 "Exponentiation of a negative base number with "
2247 "a real exponent can't be performed."));
2248 prod *=
2249 (s == d ? exponents[s] * std::pow(p[s], exponents[s] - 1) :
2250 std::pow(p[s], exponents[s]));
2251 }
2252 }
2253 r[d] = prod;
2254 }
2255
2256 return r;
2257 }
2258
2259
2260
2261 template <int dim, typename Number>
2262 void
2264 std::vector<Number> &values,
2265 const unsigned int component) const
2266 {
2267 Assert(values.size() == points.size(),
2268 ExcDimensionMismatch(values.size(), points.size()));
2269
2270 for (unsigned int i = 0; i < points.size(); ++i)
2271 values[i] = Monomial<dim, Number>::value(points[i], component);
2272 }
2273
2274
2275 template <int dim>
2276 Bessel1<dim>::Bessel1(const unsigned int order,
2277 const double wave_number,
2278 const Point<dim> center)
2279 : order(order)
2280 , wave_number(wave_number)
2281 , center(center)
2282 {
2283 Assert(wave_number >= 0., ExcMessage("wave_number must be nonnegative!"));
2284 }
2285
2286 template <int dim>
2287 double
2288 Bessel1<dim>::value(const Point<dim> &p, const unsigned int) const
2289 {
2290 Assert(dim == 2, ExcNotImplemented());
2291 const double r = p.distance(center);
2292 return std_cxx17::cyl_bessel_j(order, r * wave_number);
2293 }
2294
2295
2296 template <int dim>
2297 void
2298 Bessel1<dim>::value_list(const std::vector<Point<dim>> &points,
2299 std::vector<double> &values,
2300 const unsigned int) const
2301 {
2302 Assert(dim == 2, ExcNotImplemented());
2303 AssertDimension(points.size(), values.size());
2304 for (unsigned int k = 0; k < points.size(); ++k)
2305 {
2306 const double r = points[k].distance(center);
2307 values[k] = std_cxx17::cyl_bessel_j(order, r * wave_number);
2308 }
2309 }
2310
2311
2312 template <int dim>
2314 Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
2315 {
2316 Assert(dim == 2, ExcNotImplemented());
2317 const double r = p.distance(center);
2318 const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r;
2319 const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r;
2320
2321 const double dJn =
2322 (order == 0) ?
2323 (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2324 (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2325 std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2326 Tensor<1, dim> result;
2327 result[0] = wave_number * co * dJn;
2328 result[1] = wave_number * si * dJn;
2329 return result;
2330 }
2331
2332
2333
2334 template <int dim>
2335 void
2336 Bessel1<dim>::gradient_list(const std::vector<Point<dim>> &points,
2337 std::vector<Tensor<1, dim>> &gradients,
2338 const unsigned int) const
2339 {
2340 Assert(dim == 2, ExcNotImplemented());
2341 AssertDimension(points.size(), gradients.size());
2342 for (unsigned int k = 0; k < points.size(); ++k)
2343 {
2344 const Point<dim> &p = points[k];
2345 const double r = p.distance(center);
2346 const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r;
2347 const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r;
2348
2349 const double dJn =
2350 (order == 0) ?
2351 (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2352 (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2353 std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2354 Tensor<1, dim> &result = gradients[k];
2355 result[0] = wave_number * co * dJn;
2356 result[1] = wave_number * si * dJn;
2357 }
2358 }
2359
2360
2361
2362 namespace
2363 {
2364 // interpolate a data value from a table where ix denotes
2365 // the (lower) left endpoint of the interval to interpolate
2366 // in, and p_unit denotes the point in unit coordinates to do so.
2367 double
2368 interpolate(const Table<1, double> &data_values,
2369 const TableIndices<1> &ix,
2370 const Point<1> &xi)
2371 {
2372 return ((1 - xi[0]) * data_values[ix[0]] +
2373 xi[0] * data_values[ix[0] + 1]);
2374 }
2375
2376 double
2377 interpolate(const Table<2, double> &data_values,
2378 const TableIndices<2> &ix,
2379 const Point<2> &p_unit)
2380 {
2381 return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
2382 p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
2383 (1 - p_unit[1]) +
2384 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
2385 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
2386 p_unit[1]);
2387 }
2388
2389 double
2390 interpolate(const Table<3, double> &data_values,
2391 const TableIndices<3> &ix,
2392 const Point<3> &p_unit)
2393 {
2394 return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
2395 p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
2396 (1 - p_unit[1]) +
2397 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
2398 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
2399 p_unit[1]) *
2400 (1 - p_unit[2]) +
2401 (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
2402 p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
2403 (1 - p_unit[1]) +
2404 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
2405 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
2406 p_unit[1]) *
2407 p_unit[2]);
2408 }
2409
2410
2411 // Interpolate the gradient of a data value from a table where ix
2412 // denotes the lower left endpoint of the interval to interpolate
2413 // in, p_unit denotes the point in unit coordinates, and dx
2414 // denotes the width of the interval in each dimension.
2416 gradient_interpolate(const Table<1, double> &data_values,
2417 const TableIndices<1> &ix,
2418 const Point<1> &p_unit,
2419 const Point<1> &dx)
2420 {
2421 (void)p_unit;
2422 Tensor<1, 1> grad;
2423 grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
2424 return grad;
2425 }
2426
2427
2429 gradient_interpolate(const Table<2, double> &data_values,
2430 const TableIndices<2> &ix,
2431 const Point<2> &p_unit,
2432 const Point<2> &dx)
2433 {
2434 Tensor<1, 2> grad;
2435 double u00 = data_values[ix[0]][ix[1]],
2436 u01 = data_values[ix[0] + 1][ix[1]],
2437 u10 = data_values[ix[0]][ix[1] + 1],
2438 u11 = data_values[ix[0] + 1][ix[1] + 1];
2439
2440 grad[0] =
2441 ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
2442 grad[1] =
2443 ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
2444 return grad;
2445 }
2446
2447
2449 gradient_interpolate(const Table<3, double> &data_values,
2450 const TableIndices<3> &ix,
2451 const Point<3> &p_unit,
2452 const Point<3> &dx)
2453 {
2454 Tensor<1, 3> grad;
2455 double u000 = data_values[ix[0]][ix[1]][ix[2]],
2456 u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
2457 u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
2458 u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
2459 u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
2460 u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
2461 u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
2462 u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
2463
2464 grad[0] =
2465 ((1 - p_unit[2]) *
2466 ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
2467 p_unit[2] *
2468 ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
2469 dx[0];
2470 grad[1] =
2471 ((1 - p_unit[2]) *
2472 ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
2473 p_unit[2] *
2474 ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
2475 dx[1];
2476 grad[2] =
2477 ((1 - p_unit[1]) *
2478 ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
2479 p_unit[1] *
2480 ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
2481 dx[2];
2482
2483 return grad;
2484 }
2485 } // namespace
2486
2487
2488
2489 template <int dim>
2491 const std::array<std::vector<double>, dim> &coordinate_values,
2492 const Table<dim, double> &data_values)
2493 : coordinate_values(coordinate_values)
2494 , data_values(data_values)
2495 {
2496 for (unsigned int d = 0; d < dim; ++d)
2497 {
2498 Assert(
2499 coordinate_values[d].size() >= 2,
2500 ExcMessage(
2501 "Coordinate arrays must have at least two coordinate values!"));
2502 for (unsigned int i = 0; i < coordinate_values[d].size() - 1; ++i)
2503 Assert(
2504 coordinate_values[d][i] < coordinate_values[d][i + 1],
2505 ExcMessage(
2506 "Coordinate arrays must be sorted in strictly ascending order."));
2507
2508 Assert(data_values.size()[d] == coordinate_values[d].size(),
2509 ExcMessage(
2510 "Data and coordinate tables do not have the same size."));
2511 }
2512 }
2513
2514
2515
2516 template <int dim>
2518 std::array<std::vector<double>, dim> &&coordinate_values,
2519 Table<dim, double> &&data_values)
2520 : coordinate_values(std::move(coordinate_values))
2521 , data_values(std::move(data_values))
2522 {
2523 for (unsigned int d = 0; d < dim; ++d)
2524 {
2525 Assert(
2526 this->coordinate_values[d].size() >= 2,
2527 ExcMessage(
2528 "Coordinate arrays must have at least two coordinate values!"));
2529 for (unsigned int i = 0; i < this->coordinate_values[d].size() - 1; ++i)
2530 Assert(
2531 this->coordinate_values[d][i] < this->coordinate_values[d][i + 1],
2532 ExcMessage(
2533 "Coordinate arrays must be sorted in strictly ascending order."));
2534
2535 Assert(this->data_values.size()[d] == this->coordinate_values[d].size(),
2536 ExcMessage(
2537 "Data and coordinate tables do not have the same size."));
2538 }
2539 }
2540
2541
2542
2543 template <int dim>
2546 const Point<dim> &p) const
2547 {
2548 // find out where this data point lies, relative to the given
2549 // points. if we run all the way to the end of the range,
2550 // set the indices so that we will simply query the last of the
2551 // intervals, starting at x.size()-2 and going to x.size()-1.
2553 for (unsigned int d = 0; d < dim; ++d)
2554 {
2555 // get the index of the first element of the coordinate arrays that is
2556 // larger than p[d]
2557 ix[d] = (std::lower_bound(coordinate_values[d].begin(),
2558 coordinate_values[d].end(),
2559 p[d]) -
2560 coordinate_values[d].begin());
2561
2562 // the one we want is the index of the coordinate to the left, however,
2563 // so decrease it by one (unless we have a point to the left of all, in
2564 // which case we stay where we are; the formulas below are made in a way
2565 // that allow us to extend the function by a constant value)
2566 //
2567 // to make this work, if we got coordinate_values[d].end(), we actually
2568 // have to consider the last box which has index size()-2
2569 if (ix[d] == coordinate_values[d].size())
2570 ix[d] = coordinate_values[d].size() - 2;
2571 else if (ix[d] > 0)
2572 --ix[d];
2573 }
2574
2575 return ix;
2576 }
2577
2578
2579
2580 template <int dim>
2581 std::size_t
2583 {
2584 return sizeof(*this) +
2585 MemoryConsumption::memory_consumption(coordinate_values) -
2586 sizeof(coordinate_values) +
2588 sizeof(data_values);
2589 }
2590
2591
2592
2593 template <int dim>
2594 const Table<dim, double> &
2596 {
2597 return data_values;
2598 }
2599
2600
2601
2602 template <int dim>
2603 double
2605 const Point<dim> &p,
2606 const unsigned int component) const
2607 {
2608 (void)component;
2609 Assert(
2610 component == 0,
2611 ExcMessage(
2612 "This is a scalar function object, the component can only be zero."));
2613
2614 // find the index in the data table of the cell containing the input point
2615 const TableIndices<dim> ix = table_index_of_point(p);
2616
2617 // now compute the relative point within the interval/rectangle/box
2618 // defined by the point coordinates found above. truncate below and
2619 // above to accommodate points that may lie outside the range
2620 Point<dim> p_unit;
2621 for (unsigned int d = 0; d < dim; ++d)
2622 p_unit[d] = std::max(std::min((p[d] - coordinate_values[d][ix[d]]) /
2623 (coordinate_values[d][ix[d] + 1] -
2624 coordinate_values[d][ix[d]]),
2625 1.),
2626 0.);
2627
2628 return interpolate(data_values, ix, p_unit);
2629 }
2630
2631
2632
2633 template <int dim>
2636 const Point<dim> &p,
2637 const unsigned int component) const
2638 {
2639 (void)component;
2640 Assert(
2641 component == 0,
2642 ExcMessage(
2643 "This is a scalar function object, the component can only be zero."));
2644
2645 // find out where this data point lies
2646 const TableIndices<dim> ix = table_index_of_point(p);
2647
2648 Point<dim> dx;
2649 for (unsigned int d = 0; d < dim; ++d)
2650 dx[d] = coordinate_values[d][ix[d] + 1] - coordinate_values[d][ix[d]];
2651
2652 Point<dim> p_unit;
2653 for (unsigned int d = 0; d < dim; ++d)
2654 p_unit[d] =
2655 std::max(std::min((p[d] - coordinate_values[d][ix[d]]) / dx[d], 1.),
2656 0.0);
2657
2658 return gradient_interpolate(data_values, ix, p_unit, dx);
2659 }
2660
2661
2662
2663 template <int dim>
2665 const std::array<std::pair<double, double>, dim> &interval_endpoints,
2666 const std::array<unsigned int, dim> &n_subintervals,
2667 const Table<dim, double> &data_values)
2668 : interval_endpoints(interval_endpoints)
2669 , n_subintervals(n_subintervals)
2670 , data_values(data_values)
2671 {
2672 for (unsigned int d = 0; d < dim; ++d)
2673 {
2674 Assert(n_subintervals[d] >= 1,
2675 ExcMessage("There needs to be at least one subinterval in each "
2676 "coordinate direction."));
2678 ExcMessage("The interval in each coordinate direction needs "
2679 "to have positive size"));
2680 Assert(data_values.size()[d] == n_subintervals[d] + 1,
2681 ExcMessage("The data table does not have the correct size."));
2682 }
2683 }
2684
2685
2686
2687 template <int dim>
2689 std::array<std::pair<double, double>, dim> &&interval_endpoints,
2690 std::array<unsigned int, dim> &&n_subintervals,
2691 Table<dim, double> &&data_values)
2692 : interval_endpoints(std::move(interval_endpoints))
2693 , n_subintervals(std::move(n_subintervals))
2694 , data_values(std::move(data_values))
2695 {
2696 for (unsigned int d = 0; d < dim; ++d)
2697 {
2698 Assert(this->n_subintervals[d] >= 1,
2699 ExcMessage("There needs to be at least one subinterval in each "
2700 "coordinate direction."));
2702 this->interval_endpoints[d].second,
2703 ExcMessage("The interval in each coordinate direction needs "
2704 "to have positive size"));
2705 Assert(this->data_values.size()[d] == this->n_subintervals[d] + 1,
2706 ExcMessage("The data table does not have the correct size."));
2707 }
2708 }
2709
2710
2711
2712 template <int dim>
2713 double
2715 const unsigned int component) const
2716 {
2717 (void)component;
2718 Assert(
2719 component == 0,
2720 ExcMessage(
2721 "This is a scalar function object, the component can only be zero."));
2722
2723 // find out where this data point lies, relative to the given
2724 // subdivision points
2726 for (unsigned int d = 0; d < dim; ++d)
2727 {
2728 const double delta_x =
2729 ((interval_endpoints[d].second - interval_endpoints[d].first) /
2730 n_subintervals[d]);
2731 if (p[d] <= interval_endpoints[d].first)
2732 ix[d] = 0;
2733 else if (p[d] >= interval_endpoints[d].second - delta_x)
2734 ix[d] = n_subintervals[d] - 1;
2735 else
2736 ix[d] = static_cast<unsigned int>(
2737 (p[d] - interval_endpoints[d].first) / delta_x);
2738 }
2739
2740 // now compute the relative point within the interval/rectangle/box
2741 // defined by the point coordinates found above. truncate below and
2742 // above to accommodate points that may lie outside the range
2743 Point<dim> p_unit;
2744 for (unsigned int d = 0; d < dim; ++d)
2745 {
2746 const double delta_x =
2747 ((interval_endpoints[d].second - interval_endpoints[d].first) /
2748 n_subintervals[d]);
2749 p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
2750 ix[d] * delta_x) /
2751 delta_x,
2752 1.),
2753 0.);
2754 }
2755
2756 return interpolate(data_values, ix, p_unit);
2757 }
2758
2759
2760
2761 template <int dim>
2764 const unsigned int component) const
2765 {
2766 (void)component;
2767 Assert(
2768 component == 0,
2769 ExcMessage(
2770 "This is a scalar function object, the component can only be zero."));
2771
2772 // find out where this data point lies, relative to the given
2773 // subdivision points
2775 for (unsigned int d = 0; d < dim; ++d)
2776 {
2777 const double delta_x = ((this->interval_endpoints[d].second -
2778 this->interval_endpoints[d].first) /
2779 this->n_subintervals[d]);
2780 if (p[d] <= this->interval_endpoints[d].first)
2781 ix[d] = 0;
2782 else if (p[d] >= this->interval_endpoints[d].second - delta_x)
2783 ix[d] = this->n_subintervals[d] - 1;
2784 else
2785 ix[d] = static_cast<unsigned int>(
2786 (p[d] - this->interval_endpoints[d].first) / delta_x);
2787 }
2788
2789 // now compute the relative point within the interval/rectangle/box
2790 // defined by the point coordinates found above. truncate below and
2791 // above to accommodate points that may lie outside the range
2792 Point<dim> p_unit;
2793 Point<dim> delta_x;
2794 for (unsigned int d = 0; d < dim; ++d)
2795 {
2796 delta_x[d] = ((this->interval_endpoints[d].second -
2797 this->interval_endpoints[d].first) /
2798 this->n_subintervals[d]);
2799 p_unit[d] =
2800 std::max(std::min((p[d] - this->interval_endpoints[d].first -
2801 ix[d] * delta_x[d]) /
2802 delta_x[d],
2803 1.),
2804 0.);
2805 }
2806
2807 return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
2808 }
2809
2810
2811
2812 template <int dim>
2813 std::size_t
2815 {
2816 return sizeof(*this) + data_values.memory_consumption() -
2817 sizeof(data_values);
2818 }
2819
2820
2821
2822 template <int dim>
2823 const Table<dim, double> &
2825 {
2826 return data_values;
2827 }
2828
2829
2830
2831 /* ---------------------- Polynomial ----------------------- */
2832
2833
2834
2835 template <int dim>
2837 const std::vector<double> &coefficients)
2838 : Function<dim>(1)
2839 , exponents(exponents)
2840 , coefficients(coefficients)
2841 {
2842 Assert(exponents.n_rows() == coefficients.size(),
2843 ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
2844 Assert(exponents.n_cols() == dim,
2845 ExcDimensionMismatch(exponents.n_cols(), dim));
2846 }
2847
2848
2849
2850 template <int dim>
2851 double
2853 const unsigned int component) const
2854 {
2855 (void)component;
2856 AssertIndexRange(component, 1);
2857
2858 double sum = 0;
2859 for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2860 {
2861 double prod = 1;
2862 for (unsigned int s = 0; s < dim; ++s)
2863 {
2864 if (p[s] < 0)
2865 Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
2866 ExcMessage("Exponentiation of a negative base number with "
2867 "a real exponent can't be performed."));
2868 prod *= std::pow(p[s], exponents[monom][s]);
2869 }
2870 sum += coefficients[monom] * prod;
2871 }
2872 return sum;
2873 }
2874
2875
2876
2877 template <int dim>
2878 void
2879 Polynomial<dim>::value_list(const std::vector<Point<dim>> &points,
2880 std::vector<double> &values,
2881 const unsigned int component) const
2882 {
2883 Assert(values.size() == points.size(),
2884 ExcDimensionMismatch(values.size(), points.size()));
2885
2886 for (unsigned int i = 0; i < points.size(); ++i)
2887 values[i] = Polynomial<dim>::value(points[i], component);
2888 }
2889
2890
2891
2892 template <int dim>
2895 const unsigned int component) const
2896 {
2897 (void)component;
2898 AssertIndexRange(component, 1);
2899
2901
2902 for (unsigned int d = 0; d < dim; ++d)
2903 {
2904 double sum = 0;
2905
2906 for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2907 {
2908 double prod = 1;
2909 for (unsigned int s = 0; s < dim; ++s)
2910 {
2911 if ((s == d) && (exponents[monom][s] == 0) && (p[s] == 0))
2912 {
2913 prod = 0;
2914 break;
2915 }
2916 else
2917 {
2918 if (p[s] < 0)
2919 Assert(std::floor(exponents[monom][s]) ==
2920 exponents[monom][s],
2921 ExcMessage(
2922 "Exponentiation of a negative base number with "
2923 "a real exponent can't be performed."));
2924 prod *=
2925 (s == d ? exponents[monom][s] *
2926 std::pow(p[s], exponents[monom][s] - 1) :
2927 std::pow(p[s], exponents[monom][s]));
2928 }
2929 }
2930 sum += coefficients[monom] * prod;
2931 }
2932 r[d] = sum;
2933 }
2934 return r;
2935 }
2936
2937
2938
2939 template <int dim>
2940 std::size_t
2942 {
2943 return sizeof(*this) + exponents.memory_consumption() - sizeof(exponents) +
2945 sizeof(coefficients);
2946 }
2947
2948 template <int dim>
2950 : Function<dim>(dim)
2951 , T(T)
2952 {
2953 AssertThrow(dim > 1, ExcNotImplemented());
2954 }
2955
2956
2957 template <int dim>
2958 void
2960 Vector<double> &values) const
2961 {
2962 const double pi_x = numbers::PI * point[0];
2963 const double pi_y = numbers::PI * point[1];
2964 const double pi_t = numbers::PI / T * this->get_time();
2965
2966 values[0] = -2 * std::cos(pi_t) *
2967 Utilities::fixed_power<2>(std::sin(pi_x)) * std::sin(pi_y) *
2968 std::cos(pi_y);
2969 values[1] = +2 * std::cos(pi_t) *
2970 Utilities::fixed_power<2>(std::sin(pi_y)) * std::sin(pi_x) *
2971 std::cos(pi_x);
2972
2973 if (dim == 3)
2974 values[2] = 0;
2975 }
2976
2977
2978 // explicit instantiations
2979 template class SquareFunction<1>;
2980 template class SquareFunction<2>;
2981 template class SquareFunction<3>;
2982 template class Q1WedgeFunction<1>;
2983 template class Q1WedgeFunction<2>;
2984 template class Q1WedgeFunction<3>;
2985 template class PillowFunction<1>;
2986 template class PillowFunction<2>;
2987 template class PillowFunction<3>;
2988 template class CosineFunction<1>;
2989 template class CosineFunction<2>;
2990 template class CosineFunction<3>;
2991 template class CosineGradFunction<1>;
2992 template class CosineGradFunction<2>;
2993 template class CosineGradFunction<3>;
2994 template class ExpFunction<1>;
2995 template class ExpFunction<2>;
2996 template class ExpFunction<3>;
2997 template class JumpFunction<1>;
2998 template class JumpFunction<2>;
2999 template class JumpFunction<3>;
3000 template class FourierCosineFunction<1>;
3001 template class FourierCosineFunction<2>;
3002 template class FourierCosineFunction<3>;
3003 template class FourierSineFunction<1>;
3004 template class FourierSineFunction<2>;
3005 template class FourierSineFunction<3>;
3006 template class FourierCosineSum<1>;
3007 template class FourierCosineSum<2>;
3008 template class FourierCosineSum<3>;
3009 template class FourierSineSum<1>;
3010 template class FourierSineSum<2>;
3011 template class FourierSineSum<3>;
3012 template class SlitSingularityFunction<2>;
3013 template class SlitSingularityFunction<3>;
3014 template class Monomial<1>;
3015 template class Monomial<2>;
3016 template class Monomial<3>;
3017 template class Monomial<1, float>;
3018 template class Monomial<2, float>;
3019 template class Monomial<3, float>;
3020 template class Bessel1<1>;
3021 template class Bessel1<2>;
3022 template class Bessel1<3>;
3026 template class InterpolatedUniformGridData<1>;
3027 template class InterpolatedUniformGridData<2>;
3028 template class InterpolatedUniformGridData<3>;
3029 template class Polynomial<1>;
3030 template class Polynomial<2>;
3031 template class Polynomial<3>;
3032 template class RayleighKotheVortex<1>;
3033 template class RayleighKotheVortex<2>;
3034 template class RayleighKotheVortex<3>;
3035} // namespace Functions
3036
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double value(const Point< dim > &points, const unsigned int component=0) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Bessel1(const unsigned int order, const double wave_number, const Point< dim > center=Point< dim >())
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
CosineFunction(const unsigned int n_components=1)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
virtual void hessian_list(const std::vector< Point< dim > > &points, std::vector< SymmetricTensor< 2, dim > > &hessians, const unsigned int component=0) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component) const override
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component) const override
virtual double value(const Point< dim > &p, const unsigned int component) const override
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component) const override
virtual void vector_gradient_list(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
FourierCosineFunction(const Tensor< 1, dim > &fourier_coefficients)
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
const std::vector< Point< dim > > fourier_coefficients
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
FourierCosineSum(const std::vector< Point< dim > > &fourier_coefficients, const std::vector< double > &weights)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
const std::vector< double > weights
FourierSineFunction(const Tensor< 1, dim > &fourier_coefficients)
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
const std::vector< Point< dim > > fourier_coefficients
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
FourierSineSum(const std::vector< Point< dim > > &fourier_coefficients, const std::vector< double > &weights)
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
const std::vector< double > weights
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
InterpolatedTensorProductGridData(const std::array< std::vector< double >, dim > &coordinate_values, const Table< dim, double > &data_values)
const Table< dim, double > & get_data() const
TableIndices< dim > table_index_of_point(const Point< dim > &p) const
virtual std::size_t memory_consumption() const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
const std::array< std::vector< double >, dim > coordinate_values
const Table< dim, double > & get_data() const
const std::array< unsigned int, dim > n_subintervals
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
InterpolatedUniformGridData(const std::array< std::pair< double, double >, dim > &interval_endpoints, const std::array< unsigned int, dim > &n_subintervals, const Table< dim, double > &data_values)
const Table< dim, double > data_values
const std::array< std::pair< double, double >, dim > interval_endpoints
virtual std::size_t memory_consumption() const override
const Point< dim > direction
virtual std::size_t memory_consumption() const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
JumpFunction(const Point< dim > &direction, const double steepness)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component=0) const override
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
virtual double value(const Point< 2 > &p, const unsigned int component=0) const override
virtual void vector_gradient_list(const std::vector< Point< 2 > > &, std::vector< std::vector< Tensor< 1, 2 > > > &) const override
virtual void gradient_list(const std::vector< Point< 2 > > &points, std::vector< Tensor< 1, 2 > > &gradients, const unsigned int component=0) const override
virtual void vector_value_list(const std::vector< Point< 2 > > &points, std::vector< Vector< double > > &values) const override
virtual void vector_gradient_list(const std::vector< Point< 2 > > &, std::vector< std::vector< Tensor< 1, 2 > > > &) const override
virtual double laplacian(const Point< 2 > &p, const unsigned int component) const override
virtual void value_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component) const override
virtual void vector_value_list(const std::vector< Point< 2 > > &points, std::vector< Vector< double > > &values) const override
virtual void laplacian_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component) const override
virtual void gradient_list(const std::vector< Point< 2 > > &points, std::vector< Tensor< 1, 2 > > &gradients, const unsigned int component) const override
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component) const override
virtual double value(const Point< 2 > &p, const unsigned int component) const override
virtual Number value(const Point< dim > &p, const unsigned int component=0) const override
virtual Tensor< 1, dim, Number > gradient(const Point< dim > &p, const unsigned int component=0) const override
Monomial(const Tensor< 1, dim, Number > &exponents, const unsigned int n_components=1)
virtual void vector_value(const Point< dim > &p, Vector< Number > &values) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< Number > &values, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
PillowFunction(const double offset=0.)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
const Table< 2, double > exponents
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual std::size_t memory_consumption() const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
const std::vector< double > coefficients
Polynomial(const Table< 2, double > &exponents, const std::vector< double > &coefficients)
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_gradient_list(const std::vector< Point< dim > > &, std::vector< std::vector< Tensor< 1, dim > > > &) const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
RayleighKotheVortex(const double T=1.0)
virtual void vector_value(const Point< dim > &point, Vector< double > &values) const override
virtual void gradient_list(const std::vector< Point< 2 > > &points, std::vector< Tensor< 1, 2 > > &gradients, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void vector_gradient_list(const std::vector< Point< 2 > > &, std::vector< std::vector< Tensor< 1, 2 > > > &) const override
virtual double value(const Point< 2 > &p, const unsigned int component=0) const override
virtual void vector_value_list(const std::vector< Point< 2 > > &points, std::vector< Vector< double > > &values) const override
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void vector_gradient_list(const std::vector< Point< dim > > &, std::vector< std::vector< Tensor< 1, dim > > > &) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim > > &gradient) const override
Definition point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr numbers::NumberTraits< Number >::real_type square() const
static constexpr std::size_t memory_consumption()
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > center
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
static constexpr double PI_2
Definition numbers.h:264
static constexpr double PI
Definition numbers.h:259
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)