deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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
15#include <deal.II/base/function_bessel.h>
16#include <deal.II/base/function_lib.h>
17#include <deal.II/base/numbers.h>
18#include <deal.II/base/point.h>
19#include <deal.II/base/std_cxx17/cmath.h>
20#include <deal.II/base/tensor.h>
21
22#include <deal.II/lac/vector.h>
23
24#include <cmath>
25
26DEAL_II_NAMESPACE_OPEN
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
41 SquareFunction<dim>::vector_value(const Point<dim> &p,
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
76 SquareFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
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>
90 Tensor<1, 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
99 SquareFunction<dim>::vector_gradient(
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
111 SquareFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
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
138 Q1WedgeFunction<dim>::value_list(const std::vector<Point<dim>> &points,
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
156 Q1WedgeFunction<dim>::vector_value_list(
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
184 Q1WedgeFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
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>
199 Tensor<1, 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
213 Q1WedgeFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
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
231 Q1WedgeFunction<dim>::vector_gradient_list(
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>
253 PillowFunction<dim>::PillowFunction(const double offset)
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:
272 DEAL_II_NOT_IMPLEMENTED();
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:
303 DEAL_II_NOT_IMPLEMENTED();
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:
325 DEAL_II_NOT_IMPLEMENTED();
326 }
327 return 0.;
328 }
329
330 template <int dim>
331 void
332 PillowFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
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:
356 DEAL_II_NOT_IMPLEMENTED();
357 }
358 }
359 }
360
361 template <int dim>
362 Tensor<1, 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:
381 DEAL_II_NOT_IMPLEMENTED();
382 }
383 return result;
384 }
385
386 template <int dim>
387 void
388 PillowFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
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:
416 DEAL_II_NOT_IMPLEMENTED();
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:
446 DEAL_II_NOT_IMPLEMENTED();
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
467 CosineFunction<dim>::vector_value_list(
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:
502 DEAL_II_NOT_IMPLEMENTED();
503 }
504 return 0.;
505 }
506
507 template <int dim>
508 void
509 CosineFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
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>
521 Tensor<1, 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:
548 DEAL_II_NOT_IMPLEMENTED();
549 }
550 return result;
551 }
552
553 template <int dim>
554 void
555 CosineFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
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] =
580 -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) *
581 std::cos(numbers::PI_2 * p[1]) * std::cos(numbers::PI_2 * p[2]);
582 gradients[i][1] =
583 -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
584 std::sin(numbers::PI_2 * p[1]) * std::cos(numbers::PI_2 * p[2]);
585 gradients[i][2] =
586 -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) *
587 std::cos(numbers::PI_2 * p[1]) * std::sin(numbers::PI_2 * p[2]);
588 break;
589 default:
590 DEAL_II_NOT_IMPLEMENTED();
591 }
592 }
593 }
594
595 template <int dim>
596 SymmetricTensor<2, 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
601 SymmetricTensor<2, dim> result;
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:
644 DEAL_II_NOT_IMPLEMENTED();
645 }
646 return result;
647 }
648
649 template <int dim>
650 void
651 CosineFunction<dim>::hessian_list(
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:
706 DEAL_II_NOT_IMPLEMENTED();
707 }
708 }
709 }
710
711 //--------------------------------------------------------------------
712
713 template <int dim>
717
718
719 template <int dim>
720 double
721 CosineGradFunction<dim>::value(const Point<dim> &p,
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:
739 DEAL_II_NOT_IMPLEMENTED();
740 }
741 return 0.;
742 }
743
744
745 template <int dim>
746 void
747 CosineGradFunction<dim>::vector_value(const Point<dim> &p,
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:
774 DEAL_II_NOT_IMPLEMENTED();
775 }
776 }
777
778
779 template <int dim>
780 void
781 CosineGradFunction<dim>::value_list(const std::vector<Point<dim>> &points,
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:
809 DEAL_II_NOT_IMPLEMENTED();
810 }
811 }
812 }
813
814
815 template <int dim>
816 void
817 CosineGradFunction<dim>::vector_value_list(
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:
850 DEAL_II_NOT_IMPLEMENTED();
851 }
852 }
853 }
854
855
856 template <int dim>
857 double
858 CosineGradFunction<dim>::laplacian(const Point<dim> &p,
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>
866 Tensor<1, dim>
867 CosineGradFunction<dim>::gradient(const Point<dim> &p,
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:
899 DEAL_II_NOT_IMPLEMENTED();
900 }
901 return result;
902 }
903
904
905 template <int dim>
906 void
907 CosineGradFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
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:
946 DEAL_II_NOT_IMPLEMENTED();
947 }
948 }
949 }
950
951
952 template <int dim>
953 void
954 CosineGradFunction<dim>::vector_gradient_list(
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:
1008 DEAL_II_NOT_IMPLEMENTED();
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:
1029 DEAL_II_NOT_IMPLEMENTED();
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:
1058 DEAL_II_NOT_IMPLEMENTED();
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:
1076 DEAL_II_NOT_IMPLEMENTED();
1077 }
1078 return 0.;
1079 }
1080
1081 template <int dim>
1082 void
1083 ExpFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
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:
1105 DEAL_II_NOT_IMPLEMENTED();
1106 }
1107 }
1108 }
1109
1110 template <int dim>
1111 Tensor<1, 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:
1130 DEAL_II_NOT_IMPLEMENTED();
1131 }
1132 return result;
1133 }
1134
1135 template <int dim>
1136 void
1137 ExpFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
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:
1163 DEAL_II_NOT_IMPLEMENTED();
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
1216 LSingularityFunction::vector_value_list(
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
1254 LSingularityFunction::laplacian_list(const std::vector<Point<2>> &points,
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
1267 Tensor<1, 2>
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
1288 LSingularityFunction::gradient_list(const std::vector<Point<2>> &points,
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
1315 LSingularityFunction::vector_gradient_list(
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
1368 LSingularityGradFunction::value_list(const std::vector<Point<2>> &points,
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
1393 LSingularityGradFunction::vector_value_list(
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
1420 LSingularityGradFunction::laplacian(const Point<2> &,
1421 const unsigned int) const
1422 {
1423 return 0.;
1424 }
1425
1426
1427 void
1428 LSingularityGradFunction::laplacian_list(const std::vector<Point<2>> &points,
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
1441 Tensor<1, 2>
1442 LSingularityGradFunction::gradient(const Point<2> & /*p*/,
1443 const unsigned int /*component*/) const
1444 {
1445 DEAL_II_NOT_IMPLEMENTED();
1446 return {};
1447 }
1448
1449
1450 void
1451 LSingularityGradFunction::gradient_list(
1452 const std::vector<Point<2>> & /*points*/,
1453 std::vector<Tensor<1, 2>> & /*gradients*/,
1454 const unsigned int /*component*/) const
1455 {
1456 DEAL_II_NOT_IMPLEMENTED();
1457 }
1458
1459
1460 void
1461 LSingularityGradFunction::vector_gradient_list(
1462 const std::vector<Point<2>> & /*points*/,
1463 std::vector<std::vector<Tensor<1, 2>>> & /*gradients*/) const
1464 {
1465 DEAL_II_NOT_IMPLEMENTED();
1466 }
1467
1468 //--------------------------------------------------------------------
1469
1470 template <int dim>
1471 double
1472 SlitSingularityFunction<dim>::value(const Point<dim> &p,
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
1487 SlitSingularityFunction<dim>::value_list(
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
1510 SlitSingularityFunction<dim>::vector_value_list(
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
1535 SlitSingularityFunction<dim>::laplacian(const Point<dim> &,
1536 const unsigned int) const
1537 {
1538 return 0.;
1539 }
1540
1541
1542 template <int dim>
1543 void
1544 SlitSingularityFunction<dim>::laplacian_list(
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>
1558 Tensor<1, dim>
1559 SlitSingularityFunction<dim>::gradient(const Point<dim> &p,
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
1580 SlitSingularityFunction<dim>::gradient_list(
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
1609 SlitSingularityFunction<dim>::vector_gradient_list(
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
1642 SlitHyperSingularityFunction::value(const Point<2> &p,
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
1656 SlitHyperSingularityFunction::value_list(const std::vector<Point<2>> &points,
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
1677 SlitHyperSingularityFunction::vector_value_list(
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
1701 SlitHyperSingularityFunction::laplacian(const Point<2> &,
1702 const unsigned int) const
1703 {
1704 return 0.;
1705 }
1706
1707
1708 void
1709 SlitHyperSingularityFunction::laplacian_list(
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
1722 Tensor<1, 2>
1723 SlitHyperSingularityFunction::gradient(const Point<2> &p,
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
1744 SlitHyperSingularityFunction::gradient_list(
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
1771 SlitHyperSingularityFunction::vector_gradient_list(
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>
1801 JumpFunction<dim>::JumpFunction(const Point<dim> &direction,
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 default:
1815 angle = std::numeric_limits<double>::signaling_NaN();
1816 DEAL_II_NOT_IMPLEMENTED();
1817 }
1818 sine = std::sin(angle);
1819 cosine = std::cos(angle);
1820 }
1821
1822
1823
1824 template <int dim>
1825 double
1826 JumpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1827 {
1828 const double x = steepness * (-cosine * p[0] + sine * p[1]);
1829 return -std::atan(x);
1830 }
1831
1832
1833
1834 template <int dim>
1835 void
1836 JumpFunction<dim>::value_list(const std::vector<Point<dim>> &p,
1837 std::vector<double> &values,
1838 const unsigned int) const
1839 {
1840 Assert(values.size() == p.size(),
1841 ExcDimensionMismatch(values.size(), p.size()));
1842
1843 for (unsigned int i = 0; i < p.size(); ++i)
1844 {
1845 const double x = steepness * (-cosine * p[i][0] + sine * p[i][1]);
1846 values[i] = -std::atan(x);
1847 }
1848 }
1849
1850
1851 template <int dim>
1852 double
1853 JumpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1854 {
1855 const double x = steepness * (-cosine * p[0] + sine * p[1]);
1856 const double r = 1 + x * x;
1857 return 2 * steepness * steepness * x / (r * r);
1858 }
1859
1860
1861 template <int dim>
1862 void
1863 JumpFunction<dim>::laplacian_list(const std::vector<Point<dim>> &p,
1864 std::vector<double> &values,
1865 const unsigned int) const
1866 {
1867 Assert(values.size() == p.size(),
1868 ExcDimensionMismatch(values.size(), p.size()));
1869
1870 double f = 2 * steepness * steepness;
1871
1872 for (unsigned int i = 0; i < p.size(); ++i)
1873 {
1874 const double x = steepness * (-cosine * p[i][0] + sine * p[i][1]);
1875 const double r = 1 + x * x;
1876 values[i] = f * x / (r * r);
1877 }
1878 }
1879
1880
1881
1882 template <int dim>
1883 Tensor<1, dim>
1884 JumpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1885 {
1886 const double x = steepness * (-cosine * p[0] + sine * p[1]);
1887 const double r = -steepness * (1 + x * x);
1888 Tensor<1, dim> erg;
1889 erg[0] = cosine * r;
1890 erg[1] = sine * r;
1891 return erg;
1892 }
1893
1894
1895
1896 template <int dim>
1897 void
1898 JumpFunction<dim>::gradient_list(const std::vector<Point<dim>> &p,
1899 std::vector<Tensor<1, dim>> &gradients,
1900 const unsigned int) const
1901 {
1902 Assert(gradients.size() == p.size(),
1903 ExcDimensionMismatch(gradients.size(), p.size()));
1904
1905 for (unsigned int i = 0; i < p.size(); ++i)
1906 {
1907 const double x = steepness * (cosine * p[i][0] + sine * p[i][1]);
1908 const double r = -steepness * (1 + x * x);
1909 gradients[i][0] = cosine * r;
1910 gradients[i][1] = sine * r;
1911 }
1912 }
1913
1914
1915
1916 template <int dim>
1917 std::size_t
1918 JumpFunction<dim>::memory_consumption() const
1919 {
1920 // only simple data elements, so
1921 // use sizeof operator
1922 return sizeof(*this);
1923 }
1924
1925
1926
1927 /* ---------------------- FourierCosineFunction ----------------------- */
1928
1929
1930 template <int dim>
1931 FourierCosineFunction<dim>::FourierCosineFunction(
1932 const Tensor<1, dim> &fourier_coefficients)
1933 : Function<dim>(1)
1934 , fourier_coefficients(fourier_coefficients)
1935 {}
1936
1937
1938
1939 template <int dim>
1940 double
1941 FourierCosineFunction<dim>::value(const Point<dim> &p,
1942 const unsigned int component) const
1943 {
1944 AssertIndexRange(component, 1);
1945 return std::cos(fourier_coefficients * p);
1946 }
1947
1948
1949
1950 template <int dim>
1951 Tensor<1, dim>
1952 FourierCosineFunction<dim>::gradient(const Point<dim> &p,
1953 const unsigned int component) const
1954 {
1955 AssertIndexRange(component, 1);
1956 return -fourier_coefficients * std::sin(fourier_coefficients * p);
1957 }
1958
1959
1960
1961 template <int dim>
1962 double
1963 FourierCosineFunction<dim>::laplacian(const Point<dim> &p,
1964 const unsigned int component) const
1965 {
1966 AssertIndexRange(component, 1);
1967 return (fourier_coefficients * fourier_coefficients) *
1968 (-std::cos(fourier_coefficients * p));
1969 }
1970
1971
1972
1973 /* ---------------------- FourierSineFunction ----------------------- */
1974
1975
1976
1977 template <int dim>
1978 FourierSineFunction<dim>::FourierSineFunction(
1979 const Tensor<1, dim> &fourier_coefficients)
1980 : Function<dim>(1)
1981 , fourier_coefficients(fourier_coefficients)
1982 {}
1983
1984
1985
1986 template <int dim>
1987 double
1988 FourierSineFunction<dim>::value(const Point<dim> &p,
1989 const unsigned int component) const
1990 {
1991 AssertIndexRange(component, 1);
1992 return std::sin(fourier_coefficients * p);
1993 }
1994
1995
1996
1997 template <int dim>
1998 Tensor<1, dim>
1999 FourierSineFunction<dim>::gradient(const Point<dim> &p,
2000 const unsigned int component) const
2001 {
2002 AssertIndexRange(component, 1);
2003 return fourier_coefficients * std::cos(fourier_coefficients * p);
2004 }
2005
2006
2007
2008 template <int dim>
2009 double
2010 FourierSineFunction<dim>::laplacian(const Point<dim> &p,
2011 const unsigned int component) const
2012 {
2013 AssertIndexRange(component, 1);
2014 return (fourier_coefficients * fourier_coefficients) *
2015 (-std::sin(fourier_coefficients * p));
2016 }
2017
2018
2019
2020 /* ---------------------- FourierSineSum ----------------------- */
2021
2022
2023
2024 template <int dim>
2025 FourierSineSum<dim>::FourierSineSum(
2026 const std::vector<Point<dim>> &fourier_coefficients,
2027 const std::vector<double> &weights)
2028 : Function<dim>(1)
2029 , fourier_coefficients(fourier_coefficients)
2030 , weights(weights)
2031 {
2032 Assert(fourier_coefficients.size() > 0, ExcZero());
2033 Assert(fourier_coefficients.size() == weights.size(),
2034 ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2035 }
2036
2037
2038
2039 template <int dim>
2040 double
2041 FourierSineSum<dim>::value(const Point<dim> &p,
2042 const unsigned int component) const
2043 {
2044 AssertIndexRange(component, 1);
2045
2046 const unsigned int n = weights.size();
2047 double sum = 0;
2048 for (unsigned int s = 0; s < n; ++s)
2049 sum += weights[s] * std::sin(fourier_coefficients[s] * p);
2050
2051 return sum;
2052 }
2053
2054
2055
2056 template <int dim>
2057 Tensor<1, dim>
2058 FourierSineSum<dim>::gradient(const Point<dim> &p,
2059 const unsigned int component) const
2060 {
2061 AssertIndexRange(component, 1);
2062
2063 const unsigned int n = weights.size();
2064 Tensor<1, dim> sum;
2065 for (unsigned int s = 0; s < n; ++s)
2066 sum += fourier_coefficients[s] * std::cos(fourier_coefficients[s] * p);
2067
2068 return sum;
2069 }
2070
2071
2072
2073 template <int dim>
2074 double
2075 FourierSineSum<dim>::laplacian(const Point<dim> &p,
2076 const unsigned int component) const
2077 {
2078 AssertIndexRange(component, 1);
2079
2080 const unsigned int n = weights.size();
2081 double sum = 0;
2082 for (unsigned int s = 0; s < n; ++s)
2083 sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2084 std::sin(fourier_coefficients[s] * p);
2085
2086 return sum;
2087 }
2088
2089
2090
2091 /* ---------------------- FourierCosineSum ----------------------- */
2092
2093
2094
2095 template <int dim>
2096 FourierCosineSum<dim>::FourierCosineSum(
2097 const std::vector<Point<dim>> &fourier_coefficients,
2098 const std::vector<double> &weights)
2099 : Function<dim>(1)
2100 , fourier_coefficients(fourier_coefficients)
2101 , weights(weights)
2102 {
2103 Assert(fourier_coefficients.size() > 0, ExcZero());
2104 Assert(fourier_coefficients.size() == weights.size(),
2105 ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2106 }
2107
2108
2109
2110 template <int dim>
2111 double
2112 FourierCosineSum<dim>::value(const Point<dim> &p,
2113 const unsigned int component) const
2114 {
2115 AssertIndexRange(component, 1);
2116
2117 const unsigned int n = weights.size();
2118 double sum = 0;
2119 for (unsigned int s = 0; s < n; ++s)
2120 sum += weights[s] * std::cos(fourier_coefficients[s] * p);
2121
2122 return sum;
2123 }
2124
2125
2126
2127 template <int dim>
2128 Tensor<1, dim>
2129 FourierCosineSum<dim>::gradient(const Point<dim> &p,
2130 const unsigned int component) const
2131 {
2132 AssertIndexRange(component, 1);
2133
2134 const unsigned int n = weights.size();
2135 Tensor<1, dim> sum;
2136 for (unsigned int s = 0; s < n; ++s)
2137 sum -= fourier_coefficients[s] * std::sin(fourier_coefficients[s] * p);
2138
2139 return sum;
2140 }
2141
2142
2143
2144 template <int dim>
2145 double
2146 FourierCosineSum<dim>::laplacian(const Point<dim> &p,
2147 const unsigned int component) const
2148 {
2149 AssertIndexRange(component, 1);
2150
2151 const unsigned int n = weights.size();
2152 double sum = 0;
2153 for (unsigned int s = 0; s < n; ++s)
2154 sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2155 std::cos(fourier_coefficients[s] * p);
2156
2157 return sum;
2158 }
2159
2160
2161
2162 /* ---------------------- Monomial ----------------------- */
2163
2164
2165
2166 template <int dim, typename Number>
2167 Monomial<dim, Number>::Monomial(const Tensor<1, dim, Number> &exponents,
2168 const unsigned int n_components)
2169 : Function<dim, Number>(n_components)
2170 , exponents(exponents)
2171 {}
2172
2173
2174
2175 template <int dim, typename Number>
2176 Number
2177 Monomial<dim, Number>::value(const Point<dim> &p,
2178 const unsigned int component) const
2179 {
2180 AssertIndexRange(component, this->n_components);
2181
2182 Number prod = 1;
2183 for (unsigned int s = 0; s < dim; ++s)
2184 {
2185 if (p[s] < 0)
2186 Assert(std::floor(exponents[s]) == exponents[s],
2187 ExcMessage("Exponentiation of a negative base number with "
2188 "a real exponent can't be performed."));
2189 prod *= std::pow(p[s], exponents[s]);
2190 }
2191 return prod;
2192 }
2193
2194
2195
2196 template <int dim, typename Number>
2197 void
2198 Monomial<dim, Number>::vector_value(const Point<dim> &p,
2199 Vector<Number> &values) const
2200 {
2201 Assert(values.size() == this->n_components,
2202 ExcDimensionMismatch(values.size(), this->n_components));
2203
2204 for (unsigned int i = 0; i < values.size(); ++i)
2205 values(i) = Monomial<dim, Number>::value(p, i);
2206 }
2207
2208
2209
2210 template <int dim, typename Number>
2211 Tensor<1, dim, Number>
2212 Monomial<dim, Number>::gradient(const Point<dim> &p,
2213 const unsigned int component) const
2214 {
2215 AssertIndexRange(component, 1);
2216
2217 Tensor<1, dim> r;
2218 for (unsigned int d = 0; d < dim; ++d)
2219 {
2220 double prod = 1;
2221 for (unsigned int s = 0; s < dim; ++s)
2222 {
2223 if ((s == d) && (exponents[s] == 0) && (p[s] == 0))
2224 {
2225 prod = 0;
2226 break;
2227 }
2228 else
2229 {
2230 if (p[s] < 0)
2231 Assert(std::floor(exponents[s]) == exponents[s],
2232 ExcMessage(
2233 "Exponentiation of a negative base number with "
2234 "a real exponent can't be performed."));
2235 prod *=
2236 (s == d ? exponents[s] * std::pow(p[s], exponents[s] - 1) :
2237 std::pow(p[s], exponents[s]));
2238 }
2239 }
2240 r[d] = prod;
2241 }
2242
2243 return r;
2244 }
2245
2246
2247
2248 template <int dim, typename Number>
2249 void
2250 Monomial<dim, Number>::value_list(const std::vector<Point<dim>> &points,
2251 std::vector<Number> &values,
2252 const unsigned int component) const
2253 {
2254 Assert(values.size() == points.size(),
2255 ExcDimensionMismatch(values.size(), points.size()));
2256
2257 for (unsigned int i = 0; i < points.size(); ++i)
2258 values[i] = Monomial<dim, Number>::value(points[i], component);
2259 }
2260
2261
2262 template <int dim>
2263 Bessel1<dim>::Bessel1(const unsigned int order,
2264 const double wave_number,
2265 const Point<dim> center)
2266 : order(order)
2267 , wave_number(wave_number)
2268 , center(center)
2269 {
2270 Assert(wave_number >= 0., ExcMessage("wave_number must be nonnegative!"));
2271 }
2272
2273 template <int dim>
2274 double
2275 Bessel1<dim>::value(const Point<dim> &p, const unsigned int) const
2276 {
2277 Assert(dim == 2, ExcNotImplemented());
2278 const double r = p.distance(center);
2279 return std_cxx17::cyl_bessel_j(order, r * wave_number);
2280 }
2281
2282
2283 template <int dim>
2284 void
2285 Bessel1<dim>::value_list(const std::vector<Point<dim>> &points,
2286 std::vector<double> &values,
2287 const unsigned int) const
2288 {
2289 Assert(dim == 2, ExcNotImplemented());
2290 AssertDimension(points.size(), values.size());
2291 for (unsigned int k = 0; k < points.size(); ++k)
2292 {
2293 const double r = points[k].distance(center);
2294 values[k] = std_cxx17::cyl_bessel_j(order, r * wave_number);
2295 }
2296 }
2297
2298
2299 template <int dim>
2300 Tensor<1, dim>
2301 Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
2302 {
2303 Assert(dim == 2, ExcNotImplemented());
2304 const double r = p.distance(center);
2305 const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r;
2306 const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r;
2307
2308 const double dJn =
2309 (order == 0) ?
2310 (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2311 (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2312 std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2313 Tensor<1, dim> result;
2314 result[0] = wave_number * co * dJn;
2315 result[1] = wave_number * si * dJn;
2316 return result;
2317 }
2318
2319
2320
2321 template <int dim>
2322 void
2323 Bessel1<dim>::gradient_list(const std::vector<Point<dim>> &points,
2324 std::vector<Tensor<1, dim>> &gradients,
2325 const unsigned int) const
2326 {
2327 Assert(dim == 2, ExcNotImplemented());
2328 AssertDimension(points.size(), gradients.size());
2329 for (unsigned int k = 0; k < points.size(); ++k)
2330 {
2331 const Point<dim> &p = points[k];
2332 const double r = p.distance(center);
2333 const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r;
2334 const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r;
2335
2336 const double dJn =
2337 (order == 0) ?
2338 (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2339 (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2340 std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2341 Tensor<1, dim> &result = gradients[k];
2342 result[0] = wave_number * co * dJn;
2343 result[1] = wave_number * si * dJn;
2344 }
2345 }
2346
2347
2348
2349 namespace
2350 {
2351 // interpolate a data value from a table where ix denotes
2352 // the (lower) left endpoint of the interval to interpolate
2353 // in, and p_unit denotes the point in unit coordinates to do so.
2354 double
2355 interpolate(const Table<1, double> &data_values,
2356 const TableIndices<1> &ix,
2357 const Point<1> &xi)
2358 {
2359 return ((1 - xi[0]) * data_values[ix[0]] +
2360 xi[0] * data_values[ix[0] + 1]);
2361 }
2362
2363 double
2364 interpolate(const Table<2, double> &data_values,
2365 const TableIndices<2> &ix,
2366 const Point<2> &p_unit)
2367 {
2368 return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
2369 p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
2370 (1 - p_unit[1]) +
2371 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
2372 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
2373 p_unit[1]);
2374 }
2375
2376 double
2377 interpolate(const Table<3, double> &data_values,
2378 const TableIndices<3> &ix,
2379 const Point<3> &p_unit)
2380 {
2381 return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
2382 p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
2383 (1 - p_unit[1]) +
2384 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
2385 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
2386 p_unit[1]) *
2387 (1 - p_unit[2]) +
2388 (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
2389 p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
2390 (1 - p_unit[1]) +
2391 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
2392 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
2393 p_unit[1]) *
2394 p_unit[2]);
2395 }
2396
2397
2398 // Interpolate the gradient of a data value from a table where ix
2399 // denotes the lower left endpoint of the interval to interpolate
2400 // in, p_unit denotes the point in unit coordinates, and dx
2401 // denotes the width of the interval in each dimension.
2402 Tensor<1, 1>
2403 gradient_interpolate(const Table<1, double> &data_values,
2404 const TableIndices<1> &ix,
2405 const Point<1> & /*p_unit*/,
2406 const Point<1> &dx)
2407 {
2408 Tensor<1, 1> grad;
2409 grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
2410 return grad;
2411 }
2412
2413
2414 Tensor<1, 2>
2415 gradient_interpolate(const Table<2, double> &data_values,
2416 const TableIndices<2> &ix,
2417 const Point<2> &p_unit,
2418 const Point<2> &dx)
2419 {
2420 Tensor<1, 2> grad;
2421 double u00 = data_values[ix[0]][ix[1]],
2422 u01 = data_values[ix[0] + 1][ix[1]],
2423 u10 = data_values[ix[0]][ix[1] + 1],
2424 u11 = data_values[ix[0] + 1][ix[1] + 1];
2425
2426 grad[0] =
2427 ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
2428 grad[1] =
2429 ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
2430 return grad;
2431 }
2432
2433
2434 Tensor<1, 3>
2435 gradient_interpolate(const Table<3, double> &data_values,
2436 const TableIndices<3> &ix,
2437 const Point<3> &p_unit,
2438 const Point<3> &dx)
2439 {
2440 Tensor<1, 3> grad;
2441 double u000 = data_values[ix[0]][ix[1]][ix[2]],
2442 u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
2443 u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
2444 u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
2445 u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
2446 u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
2447 u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
2448 u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
2449
2450 grad[0] =
2451 ((1 - p_unit[2]) *
2452 ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
2453 p_unit[2] *
2454 ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
2455 dx[0];
2456 grad[1] =
2457 ((1 - p_unit[2]) *
2458 ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
2459 p_unit[2] *
2460 ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
2461 dx[1];
2462 grad[2] =
2463 ((1 - p_unit[1]) *
2464 ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
2465 p_unit[1] *
2466 ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
2467 dx[2];
2468
2469 return grad;
2470 }
2471 } // namespace
2472
2473
2474
2475 template <int dim>
2476 InterpolatedTensorProductGridData<dim>::InterpolatedTensorProductGridData(
2477 const std::array<std::vector<double>, dim> &coordinate_values,
2478 const Table<dim, double> &data_values)
2479 : coordinate_values(coordinate_values)
2480 , data_values(data_values)
2481 {
2482 for (unsigned int d = 0; d < dim; ++d)
2483 {
2484 Assert(
2485 coordinate_values[d].size() >= 2,
2486 ExcMessage(
2487 "Coordinate arrays must have at least two coordinate values!"));
2488 for (unsigned int i = 0; i < coordinate_values[d].size() - 1; ++i)
2489 Assert(
2490 coordinate_values[d][i] < coordinate_values[d][i + 1],
2491 ExcMessage(
2492 "Coordinate arrays must be sorted in strictly ascending order."));
2493
2494 Assert(data_values.size()[d] == coordinate_values[d].size(),
2495 ExcMessage(
2496 "Data and coordinate tables do not have the same size."));
2497 }
2498 }
2499
2500
2501
2502 template <int dim>
2503 InterpolatedTensorProductGridData<dim>::InterpolatedTensorProductGridData(
2504 std::array<std::vector<double>, dim> &&coordinate_values,
2505 Table<dim, double> &&data_values)
2506 : coordinate_values(std::move(coordinate_values))
2507 , data_values(std::move(data_values))
2508 {
2509 for (unsigned int d = 0; d < dim; ++d)
2510 {
2511 Assert(
2512 this->coordinate_values[d].size() >= 2,
2513 ExcMessage(
2514 "Coordinate arrays must have at least two coordinate values!"));
2515 for (unsigned int i = 0; i < this->coordinate_values[d].size() - 1; ++i)
2516 Assert(
2517 this->coordinate_values[d][i] < this->coordinate_values[d][i + 1],
2518 ExcMessage(
2519 "Coordinate arrays must be sorted in strictly ascending order."));
2520
2521 Assert(this->data_values.size()[d] == this->coordinate_values[d].size(),
2522 ExcMessage(
2523 "Data and coordinate tables do not have the same size."));
2524 }
2525 }
2526
2527
2528
2529 template <int dim>
2530 TableIndices<dim>
2531 InterpolatedTensorProductGridData<dim>::table_index_of_point(
2532 const Point<dim> &p) const
2533 {
2534 // find out where this data point lies, relative to the given
2535 // points. if we run all the way to the end of the range,
2536 // set the indices so that we will simply query the last of the
2537 // intervals, starting at x.size()-2 and going to x.size()-1.
2538 TableIndices<dim> ix;
2539 for (unsigned int d = 0; d < dim; ++d)
2540 {
2541 // get the index of the first element of the coordinate arrays that is
2542 // larger than p[d]
2543 ix[d] = (std::lower_bound(coordinate_values[d].begin(),
2544 coordinate_values[d].end(),
2545 p[d]) -
2546 coordinate_values[d].begin());
2547
2548 // the one we want is the index of the coordinate to the left, however,
2549 // so decrease it by one (unless we have a point to the left of all, in
2550 // which case we stay where we are; the formulas below are made in a way
2551 // that allow us to extend the function by a constant value)
2552 //
2553 // to make this work, if we got coordinate_values[d].end(), we actually
2554 // have to consider the last box which has index size()-2
2555 if (ix[d] == coordinate_values[d].size())
2556 ix[d] = coordinate_values[d].size() - 2;
2557 else if (ix[d] > 0)
2558 --ix[d];
2559 }
2560
2561 return ix;
2562 }
2563
2564
2565
2566 template <int dim>
2567 std::size_t
2568 InterpolatedTensorProductGridData<dim>::memory_consumption() const
2569 {
2570 return sizeof(*this) +
2571 MemoryConsumption::memory_consumption(coordinate_values) -
2572 sizeof(coordinate_values) +
2573 MemoryConsumption::memory_consumption(data_values) -
2574 sizeof(data_values);
2575 }
2576
2577
2578
2579 template <int dim>
2580 const Table<dim, double> &
2581 InterpolatedTensorProductGridData<dim>::get_data() const
2582 {
2583 return data_values;
2584 }
2585
2586
2587
2588 template <int dim>
2589 double
2590 InterpolatedTensorProductGridData<dim>::value(
2591 const Point<dim> &p,
2592 const unsigned int component) const
2593 {
2594 Assert(
2595 component == 0,
2596 ExcMessage(
2597 "This is a scalar function object, the component can only be zero."));
2598
2599 // find the index in the data table of the cell containing the input point
2600 const TableIndices<dim> ix = table_index_of_point(p);
2601
2602 // now compute the relative point within the interval/rectangle/box
2603 // defined by the point coordinates found above. truncate below and
2604 // above to accommodate points that may lie outside the range
2605 Point<dim> p_unit;
2606 for (unsigned int d = 0; d < dim; ++d)
2607 p_unit[d] = std::max(std::min((p[d] - coordinate_values[d][ix[d]]) /
2608 (coordinate_values[d][ix[d] + 1] -
2609 coordinate_values[d][ix[d]]),
2610 1.),
2611 0.);
2612
2613 return interpolate(data_values, ix, p_unit);
2614 }
2615
2616
2617
2618 template <int dim>
2619 Tensor<1, dim>
2620 InterpolatedTensorProductGridData<dim>::gradient(
2621 const Point<dim> &p,
2622 const unsigned int component) const
2623 {
2624 Assert(
2625 component == 0,
2626 ExcMessage(
2627 "This is a scalar function object, the component can only be zero."));
2628
2629 // find out where this data point lies
2630 const TableIndices<dim> ix = table_index_of_point(p);
2631
2632 Point<dim> dx;
2633 for (unsigned int d = 0; d < dim; ++d)
2634 dx[d] = coordinate_values[d][ix[d] + 1] - coordinate_values[d][ix[d]];
2635
2636 Point<dim> p_unit;
2637 for (unsigned int d = 0; d < dim; ++d)
2638 p_unit[d] =
2639 std::max(std::min((p[d] - coordinate_values[d][ix[d]]) / dx[d], 1.),
2640 0.0);
2641
2642 return gradient_interpolate(data_values, ix, p_unit, dx);
2643 }
2644
2645
2646
2647 template <int dim>
2648 InterpolatedUniformGridData<dim>::InterpolatedUniformGridData(
2649 const std::array<std::pair<double, double>, dim> &interval_endpoints,
2650 const std::array<unsigned int, dim> &n_subintervals,
2651 const Table<dim, double> &data_values)
2652 : interval_endpoints(interval_endpoints)
2653 , n_subintervals(n_subintervals)
2654 , data_values(data_values)
2655 {
2656 for (unsigned int d = 0; d < dim; ++d)
2657 {
2658 Assert(n_subintervals[d] >= 1,
2659 ExcMessage("There needs to be at least one subinterval in each "
2660 "coordinate direction."));
2661 Assert(interval_endpoints[d].first < interval_endpoints[d].second,
2662 ExcMessage("The interval in each coordinate direction needs "
2663 "to have positive size"));
2664 Assert(data_values.size()[d] == n_subintervals[d] + 1,
2665 ExcMessage("The data table does not have the correct size."));
2666 }
2667 }
2668
2669
2670
2671 template <int dim>
2672 InterpolatedUniformGridData<dim>::InterpolatedUniformGridData(
2673 std::array<std::pair<double, double>, dim> &&interval_endpoints,
2674 std::array<unsigned int, dim> &&n_subintervals,
2675 Table<dim, double> &&data_values)
2676 : interval_endpoints(std::move(interval_endpoints))
2677 , n_subintervals(std::move(n_subintervals))
2678 , data_values(std::move(data_values))
2679 {
2680 for (unsigned int d = 0; d < dim; ++d)
2681 {
2682 Assert(this->n_subintervals[d] >= 1,
2683 ExcMessage("There needs to be at least one subinterval in each "
2684 "coordinate direction."));
2685 Assert(this->interval_endpoints[d].first <
2686 this->interval_endpoints[d].second,
2687 ExcMessage("The interval in each coordinate direction needs "
2688 "to have positive size"));
2689 Assert(this->data_values.size()[d] == this->n_subintervals[d] + 1,
2690 ExcMessage("The data table does not have the correct size."));
2691 }
2692 }
2693
2694
2695
2696 template <int dim>
2697 double
2698 InterpolatedUniformGridData<dim>::value(const Point<dim> &p,
2699 const unsigned int component) const
2700 {
2701 Assert(
2702 component == 0,
2703 ExcMessage(
2704 "This is a scalar function object, the component can only be zero."));
2705
2706 // find out where this data point lies, relative to the given
2707 // subdivision points
2708 TableIndices<dim> ix;
2709 for (unsigned int d = 0; d < dim; ++d)
2710 {
2711 const double delta_x =
2712 ((interval_endpoints[d].second - interval_endpoints[d].first) /
2713 n_subintervals[d]);
2714 if (p[d] <= interval_endpoints[d].first)
2715 ix[d] = 0;
2716 else if (p[d] >= interval_endpoints[d].second - delta_x)
2717 ix[d] = n_subintervals[d] - 1;
2718 else
2719 ix[d] = static_cast<unsigned int>(
2720 (p[d] - interval_endpoints[d].first) / delta_x);
2721 }
2722
2723 // now compute the relative point within the interval/rectangle/box
2724 // defined by the point coordinates found above. truncate below and
2725 // above to accommodate points that may lie outside the range
2726 Point<dim> p_unit;
2727 for (unsigned int d = 0; d < dim; ++d)
2728 {
2729 const double delta_x =
2730 ((interval_endpoints[d].second - interval_endpoints[d].first) /
2731 n_subintervals[d]);
2732 p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
2733 ix[d] * delta_x) /
2734 delta_x,
2735 1.),
2736 0.);
2737 }
2738
2739 return interpolate(data_values, ix, p_unit);
2740 }
2741
2742
2743
2744 template <int dim>
2745 Tensor<1, dim>
2746 InterpolatedUniformGridData<dim>::gradient(const Point<dim> &p,
2747 const unsigned int component) const
2748 {
2749 Assert(
2750 component == 0,
2751 ExcMessage(
2752 "This is a scalar function object, the component can only be zero."));
2753
2754 // find out where this data point lies, relative to the given
2755 // subdivision points
2756 TableIndices<dim> ix;
2757 for (unsigned int d = 0; d < dim; ++d)
2758 {
2759 const double delta_x = ((this->interval_endpoints[d].second -
2760 this->interval_endpoints[d].first) /
2761 this->n_subintervals[d]);
2762 if (p[d] <= this->interval_endpoints[d].first)
2763 ix[d] = 0;
2764 else if (p[d] >= this->interval_endpoints[d].second - delta_x)
2765 ix[d] = this->n_subintervals[d] - 1;
2766 else
2767 ix[d] = static_cast<unsigned int>(
2768 (p[d] - this->interval_endpoints[d].first) / delta_x);
2769 }
2770
2771 // now compute the relative point within the interval/rectangle/box
2772 // defined by the point coordinates found above. truncate below and
2773 // above to accommodate points that may lie outside the range
2774 Point<dim> p_unit;
2775 Point<dim> delta_x;
2776 for (unsigned int d = 0; d < dim; ++d)
2777 {
2778 delta_x[d] = ((this->interval_endpoints[d].second -
2779 this->interval_endpoints[d].first) /
2780 this->n_subintervals[d]);
2781 p_unit[d] =
2782 std::max(std::min((p[d] - this->interval_endpoints[d].first -
2783 ix[d] * delta_x[d]) /
2784 delta_x[d],
2785 1.),
2786 0.);
2787 }
2788
2789 return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
2790 }
2791
2792
2793
2794 template <int dim>
2795 std::size_t
2796 InterpolatedUniformGridData<dim>::memory_consumption() const
2797 {
2798 return sizeof(*this) + data_values.memory_consumption() -
2799 sizeof(data_values);
2800 }
2801
2802
2803
2804 template <int dim>
2805 const Table<dim, double> &
2806 InterpolatedUniformGridData<dim>::get_data() const
2807 {
2808 return data_values;
2809 }
2810
2811
2812
2813 /* ---------------------- Polynomial ----------------------- */
2814
2815
2816
2817 template <int dim>
2818 Polynomial<dim>::Polynomial(const Table<2, double> &exponents,
2819 const std::vector<double> &coefficients)
2820 : Function<dim>(1)
2821 , exponents(exponents)
2822 , coefficients(coefficients)
2823 {
2824 Assert(exponents.n_rows() == coefficients.size(),
2825 ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
2826 Assert(exponents.n_cols() == dim,
2827 ExcDimensionMismatch(exponents.n_cols(), dim));
2828 }
2829
2830
2831
2832 template <int dim>
2833 double
2834 Polynomial<dim>::value(const Point<dim> &p,
2835 const unsigned int component) const
2836 {
2837 AssertIndexRange(component, 1);
2838
2839 double sum = 0;
2840 for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2841 {
2842 double prod = 1;
2843 for (unsigned int s = 0; s < dim; ++s)
2844 {
2845 if (p[s] < 0)
2846 Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
2847 ExcMessage("Exponentiation of a negative base number with "
2848 "a real exponent can't be performed."));
2849 prod *= std::pow(p[s], exponents[monom][s]);
2850 }
2851 sum += coefficients[monom] * prod;
2852 }
2853 return sum;
2854 }
2855
2856
2857
2858 template <int dim>
2859 void
2860 Polynomial<dim>::value_list(const std::vector<Point<dim>> &points,
2861 std::vector<double> &values,
2862 const unsigned int component) const
2863 {
2864 Assert(values.size() == points.size(),
2865 ExcDimensionMismatch(values.size(), points.size()));
2866
2867 for (unsigned int i = 0; i < points.size(); ++i)
2868 values[i] = Polynomial<dim>::value(points[i], component);
2869 }
2870
2871
2872
2873 template <int dim>
2874 Tensor<1, dim>
2875 Polynomial<dim>::gradient(const Point<dim> &p,
2876 const unsigned int component) const
2877 {
2878 AssertIndexRange(component, 1);
2879
2880 Tensor<1, dim> r;
2881
2882 for (unsigned int d = 0; d < dim; ++d)
2883 {
2884 double sum = 0;
2885
2886 for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2887 {
2888 double prod = 1;
2889 for (unsigned int s = 0; s < dim; ++s)
2890 {
2891 if ((s == d) && (exponents[monom][s] == 0) && (p[s] == 0))
2892 {
2893 prod = 0;
2894 break;
2895 }
2896 else
2897 {
2898 if (p[s] < 0)
2899 Assert(std::floor(exponents[monom][s]) ==
2900 exponents[monom][s],
2901 ExcMessage(
2902 "Exponentiation of a negative base number with "
2903 "a real exponent can't be performed."));
2904 prod *=
2905 (s == d ? exponents[monom][s] *
2906 std::pow(p[s], exponents[monom][s] - 1) :
2907 std::pow(p[s], exponents[monom][s]));
2908 }
2909 }
2910 sum += coefficients[monom] * prod;
2911 }
2912 r[d] = sum;
2913 }
2914 return r;
2915 }
2916
2917
2918
2919 template <int dim>
2920 std::size_t
2921 Polynomial<dim>::memory_consumption() const
2922 {
2923 return sizeof(*this) + exponents.memory_consumption() - sizeof(exponents) +
2924 MemoryConsumption::memory_consumption(coefficients) -
2925 sizeof(coefficients);
2926 }
2927
2928 template <int dim>
2929 RayleighKotheVortex<dim>::RayleighKotheVortex(const double T)
2930 : Function<dim>(dim)
2931 , T(T)
2932 {
2933 AssertThrow(dim > 1, ExcNotImplemented());
2934 }
2935
2936
2937 template <int dim>
2938 void
2939 RayleighKotheVortex<dim>::vector_value(const Point<dim> &point,
2940 Vector<double> &values) const
2941 {
2942 const double pi_x = numbers::PI * point[0];
2943 const double pi_y = numbers::PI * point[1];
2944 const double pi_t = numbers::PI / T * this->get_time();
2945
2946 values[0] = -2 * std::cos(pi_t) *
2947 Utilities::fixed_power<2>(std::sin(pi_x)) * std::sin(pi_y) *
2948 std::cos(pi_y);
2949 values[1] = +2 * std::cos(pi_t) *
2950 Utilities::fixed_power<2>(std::sin(pi_y)) * std::sin(pi_x) *
2951 std::cos(pi_x);
2952
2953 if (dim == 3)
2954 values[2] = 0;
2955 }
2956
2957
2958 // explicit instantiations
2959 template class SquareFunction<1>;
2960 template class SquareFunction<2>;
2961 template class SquareFunction<3>;
2962 template class Q1WedgeFunction<1>;
2963 template class Q1WedgeFunction<2>;
2964 template class Q1WedgeFunction<3>;
2965 template class PillowFunction<1>;
2966 template class PillowFunction<2>;
2967 template class PillowFunction<3>;
2968 template class CosineFunction<1>;
2969 template class CosineFunction<2>;
2970 template class CosineFunction<3>;
2971 template class CosineGradFunction<1>;
2972 template class CosineGradFunction<2>;
2973 template class CosineGradFunction<3>;
2974 template class ExpFunction<1>;
2975 template class ExpFunction<2>;
2976 template class ExpFunction<3>;
2977 template class JumpFunction<1>;
2978 template class JumpFunction<2>;
2979 template class JumpFunction<3>;
2980 template class FourierCosineFunction<1>;
2981 template class FourierCosineFunction<2>;
2982 template class FourierCosineFunction<3>;
2983 template class FourierSineFunction<1>;
2984 template class FourierSineFunction<2>;
2985 template class FourierSineFunction<3>;
2986 template class FourierCosineSum<1>;
2987 template class FourierCosineSum<2>;
2988 template class FourierCosineSum<3>;
2989 template class FourierSineSum<1>;
2990 template class FourierSineSum<2>;
2991 template class FourierSineSum<3>;
2992 template class SlitSingularityFunction<2>;
2993 template class SlitSingularityFunction<3>;
2994 template class Monomial<1>;
2995 template class Monomial<2>;
2996 template class Monomial<3>;
2997 template class Monomial<1, float>;
2998 template class Monomial<2, float>;
2999 template class Monomial<3, float>;
3000 template class Bessel1<1>;
3001 template class Bessel1<2>;
3002 template class Bessel1<3>;
3003 template class InterpolatedTensorProductGridData<1>;
3004 template class InterpolatedTensorProductGridData<2>;
3005 template class InterpolatedTensorProductGridData<3>;
3006 template class InterpolatedUniformGridData<1>;
3007 template class InterpolatedUniformGridData<2>;
3008 template class InterpolatedUniformGridData<3>;
3009 template class Polynomial<1>;
3010 template class Polynomial<2>;
3011 template class Polynomial<3>;
3012 template class RayleighKotheVortex<1>;
3013 template class RayleighKotheVortex<2>;
3014 template class RayleighKotheVortex<3>;
3015} // namespace Functions
3016
3017DEAL_II_NAMESPACE_CLOSE
Functions::Bessel1::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2301
Functions::Bessel1::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:2285
Functions::Bessel1::value
virtual double value(const Point< dim > &points, const unsigned int component=0) const override
Definition function_lib.cc:2275
Functions::Bessel1::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:2323
Functions::Bessel1::Bessel1
Bessel1(const unsigned int order, const double wave_number, const Point< dim > center=Point< dim >())
Definition function_lib.cc:2263
Functions::CosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:485
Functions::CosineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:432
Functions::CosineFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:509
Functions::CosineFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:453
Functions::CosineFunction::CosineFunction
CosineFunction(const unsigned int n_components=1)
Definition function_lib.cc:424
Functions::CosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:522
Functions::CosineFunction::hessian
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:597
Functions::CosineFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:467
Functions::CosineFunction::hessian_list
virtual void hessian_list(const std::vector< Point< dim > > &points, std::vector< SymmetricTensor< 2, dim > > &hessians, const unsigned int component=0) const override
Definition function_lib.cc:651
Functions::CosineFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:555
Functions::CosineGradFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component) const override
Definition function_lib.cc:781
Functions::CosineGradFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component) const override
Definition function_lib.cc:867
Functions::CosineGradFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:817
Functions::CosineGradFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component) const override
Definition function_lib.cc:907
Functions::CosineGradFunction::value
virtual double value(const Point< dim > &p, const unsigned int component) const override
Definition function_lib.cc:721
Functions::CosineGradFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition function_lib.cc:747
Functions::CosineGradFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component) const override
Definition function_lib.cc:858
Functions::CosineGradFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const override
Definition function_lib.cc:954
Functions::ExpFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1112
Functions::ExpFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1083
Functions::ExpFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1018
Functions::ExpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1065
Functions::ExpFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1036
Functions::ExpFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:1137
Functions::FourierCosineFunction::FourierCosineFunction
FourierCosineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition function_lib.cc:1931
Functions::FourierCosineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1941
Functions::FourierCosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1963
Functions::FourierCosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1952
Functions::FourierCosineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition function_lib.h:907
Functions::FourierCosineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2146
Functions::FourierCosineSum::FourierCosineSum
FourierCosineSum(const std::vector< Point< dim > > &fourier_coefficients, const std::vector< double > &weights)
Definition function_lib.cc:2096
Functions::FourierCosineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2129
Functions::FourierCosineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2112
Functions::FourierCosineSum::weights
const std::vector< double > weights
Definition function_lib.h:908
Functions::FourierSineFunction::FourierSineFunction
FourierSineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition function_lib.cc:1978
Functions::FourierSineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2010
Functions::FourierSineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1988
Functions::FourierSineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1999
Functions::FourierSineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition function_lib.h:853
Functions::FourierSineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2058
Functions::FourierSineSum::FourierSineSum
FourierSineSum(const std::vector< Point< dim > > &fourier_coefficients, const std::vector< double > &weights)
Definition function_lib.cc:2025
Functions::FourierSineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2041
Functions::FourierSineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2075
Functions::FourierSineSum::weights
const std::vector< double > weights
Definition function_lib.h:854
Functions::InterpolatedTensorProductGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2590
Functions::InterpolatedTensorProductGridData::InterpolatedTensorProductGridData
InterpolatedTensorProductGridData(const std::array< std::vector< double >, dim > &coordinate_values, const Table< dim, double > &data_values)
Definition function_lib.cc:2476
Functions::InterpolatedTensorProductGridData::get_data
const Table< dim, double > & get_data() const
Definition function_lib.cc:2581
Functions::InterpolatedTensorProductGridData::table_index_of_point
TableIndices< dim > table_index_of_point(const Point< dim > &p) const
Definition function_lib.cc:2531
Functions::InterpolatedTensorProductGridData::memory_consumption
virtual std::size_t memory_consumption() const override
Definition function_lib.cc:2568
Functions::InterpolatedTensorProductGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2620
Functions::InterpolatedTensorProductGridData::coordinate_values
const std::array< std::vector< double >, dim > coordinate_values
Definition function_lib.h:1546
Functions::InterpolatedUniformGridData::get_data
const Table< dim, double > & get_data() const
Definition function_lib.cc:2806
Functions::InterpolatedUniformGridData::n_subintervals
const std::array< unsigned int, dim > n_subintervals
Definition function_lib.h:1691
Functions::InterpolatedUniformGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2746
Functions::InterpolatedUniformGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2698
Functions::InterpolatedUniformGridData::InterpolatedUniformGridData
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)
Definition function_lib.cc:2648
Functions::InterpolatedUniformGridData::data_values
const Table< dim, double > data_values
Definition function_lib.h:1696
Functions::InterpolatedUniformGridData::interval_endpoints
const std::array< std::pair< double, double >, dim > interval_endpoints
Definition function_lib.h:1686
Functions::InterpolatedUniformGridData::memory_consumption
virtual std::size_t memory_consumption() const override
Definition function_lib.cc:2796
Functions::JumpFunction::direction
const Point< dim > direction
Definition function_lib.h:674
Functions::JumpFunction::memory_consumption
virtual std::size_t memory_consumption() const override
Definition function_lib.cc:1918
Functions::JumpFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1863
Functions::JumpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1853
Functions::JumpFunction::JumpFunction
JumpFunction(const Point< dim > &direction, const double steepness)
Definition function_lib.cc:1801
Functions::JumpFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1884
Functions::JumpFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1836
Functions::JumpFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:1898
Functions::JumpFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1826
Functions::LSingularityFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1254
Functions::LSingularityFunction::value_list
virtual void value_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1189
Functions::LSingularityFunction::gradient
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component=0) const override
Definition function_lib.cc:1268
Functions::LSingularityFunction::laplacian
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
Definition function_lib.cc:1245
Functions::LSingularityFunction::value
virtual double value(const Point< 2 > &p, const unsigned int component=0) const override
Definition function_lib.cc:1172
Functions::LSingularityFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< 2 > > &, std::vector< std::vector< Tensor< 1, 2 > > > &) const override
Definition function_lib.cc:1315
Functions::LSingularityFunction::gradient_list
virtual void gradient_list(const std::vector< Point< 2 > > &points, std::vector< Tensor< 1, 2 > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:1288
Functions::LSingularityFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< 2 > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:1216
Functions::LSingularityGradFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< 2 > > &, std::vector< std::vector< Tensor< 1, 2 > > > &) const override
Definition function_lib.cc:1461
Functions::LSingularityGradFunction::laplacian
virtual double laplacian(const Point< 2 > &p, const unsigned int component) const override
Definition function_lib.cc:1420
Functions::LSingularityGradFunction::value_list
virtual void value_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component) const override
Definition function_lib.cc:1368
Functions::LSingularityGradFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< 2 > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:1393
Functions::LSingularityGradFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component) const override
Definition function_lib.cc:1428
Functions::LSingularityGradFunction::gradient_list
virtual void gradient_list(const std::vector< Point< 2 > > &points, std::vector< Tensor< 1, 2 > > &gradients, const unsigned int component) const override
Definition function_lib.cc:1451
Functions::LSingularityGradFunction::gradient
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component) const override
Definition function_lib.cc:1442
Functions::LSingularityGradFunction::value
virtual double value(const Point< 2 > &p, const unsigned int component) const override
Definition function_lib.cc:1350
Functions::Monomial::value
virtual Number value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2177
Functions::Monomial::gradient
virtual Tensor< 1, dim, Number > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2212
Functions::Monomial::Monomial
Monomial(const Tensor< 1, dim, Number > &exponents, const unsigned int n_components=1)
Definition function_lib.cc:2167
Functions::Monomial::vector_value
virtual void vector_value(const Point< dim > &p, Vector< Number > &values) const override
Definition function_lib.cc:2198
Functions::Monomial::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< Number > &values, const unsigned int component=0) const override
Definition function_lib.cc:2250
Functions::PillowFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:332
Functions::PillowFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:312
Functions::PillowFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:388
Functions::PillowFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:260
Functions::PillowFunction::PillowFunction
PillowFunction(const double offset=0.)
Definition function_lib.cc:253
Functions::PillowFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:363
Functions::PillowFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:279
Functions::Polynomial::exponents
const Table< 2, double > exponents
Definition function_lib.h:1760
Functions::Polynomial::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:2860
Functions::Polynomial::memory_consumption
virtual std::size_t memory_consumption() const override
Definition function_lib.cc:2921
Functions::Polynomial::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2834
Functions::Polynomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2875
Functions::Polynomial::coefficients
const std::vector< double > coefficients
Definition function_lib.h:1765
Functions::Polynomial::Polynomial
Polynomial(const Table< 2, double > &exponents, const std::vector< double > &coefficients)
Definition function_lib.cc:2818
Functions::Q1WedgeFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:213
Functions::Q1WedgeFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:156
Functions::Q1WedgeFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:175
Functions::Q1WedgeFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< dim > > &, std::vector< std::vector< Tensor< 1, dim > > > &) const override
Definition function_lib.cc:231
Functions::Q1WedgeFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:184
Functions::Q1WedgeFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:128
Functions::Q1WedgeFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:138
Functions::Q1WedgeFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:200
Functions::RayleighKotheVortex::RayleighKotheVortex
RayleighKotheVortex(const double T=1.0)
Definition function_lib.cc:2929
Functions::RayleighKotheVortex::vector_value
virtual void vector_value(const Point< dim > &point, Vector< double > &values) const override
Definition function_lib.cc:2939
Functions::SlitHyperSingularityFunction::gradient_list
virtual void gradient_list(const std::vector< Point< 2 > > &points, std::vector< Tensor< 1, 2 > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:1744
Functions::SlitHyperSingularityFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1709
Functions::SlitHyperSingularityFunction::value_list
virtual void value_list(const std::vector< Point< 2 > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1656
Functions::SlitHyperSingularityFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< 2 > > &, std::vector< std::vector< Tensor< 1, 2 > > > &) const override
Definition function_lib.cc:1771
Functions::SlitHyperSingularityFunction::value
virtual double value(const Point< 2 > &p, const unsigned int component=0) const override
Definition function_lib.cc:1642
Functions::SlitHyperSingularityFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< 2 > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:1677
Functions::SlitHyperSingularityFunction::laplacian
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
Definition function_lib.cc:1701
Functions::SlitHyperSingularityFunction::gradient
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component=0) const override
Definition function_lib.cc:1723
Functions::SlitSingularityFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1559
Functions::SlitSingularityFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1535
Functions::SlitSingularityFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
Definition function_lib.cc:1510
Functions::SlitSingularityFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:1580
Functions::SlitSingularityFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1487
Functions::SlitSingularityFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:1544
Functions::SlitSingularityFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< dim > > &, std::vector< std::vector< Tensor< 1, dim > > > &) const override
Definition function_lib.cc:1609
Functions::SlitSingularityFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1472
Functions::SquareFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:33
Functions::SquareFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:68
Functions::SquareFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:51
Functions::SquareFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition function_lib.cc:41
Functions::SquareFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
Definition function_lib.cc:111
Functions::SquareFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:91
Functions::SquareFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition function_lib.cc:76
Functions::SquareFunction::vector_gradient
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim > > &gradient) const override
Definition function_lib.cc:99
Point
Definition point.h:111
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Point::square
constexpr numbers::NumberTraits< Number >::real_type square() const
Tensor::memory_consumption
static constexpr std::size_t memory_consumption()
Vector::size
virtual size_type size() const override
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
function_bessel.h
function_lib.h
second
Point< 2 > second
Definition grid_out.cc:4624
first
Point< 2 > first
Definition grid_out.cc:4623
StandardExceptions::ExcZero
static ::ExceptionBase & ExcZero()
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition exceptions.h:2070
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1780
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
size
std::size_t size
Definition mpi.cc:734
FETools::interpolate
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
MemoryConsumption::memory_consumption
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Definition memory_consumption.h:271
numbers::PI_2
static constexpr double PI_2
Definition numbers.h:259
numbers::PI
static constexpr double PI
Definition numbers.h:254
std
STL namespace.
std::exp
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6829
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
std::cos
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6609
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6589
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition vectorization.h:6885
std::atan
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
Definition vectorization.h:6689
Generated by doxygen 1.9.8