Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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 (void)component;
1945 AssertIndexRange(component, 1);
1946 return std::cos(fourier_coefficients * p);
1947 }
1948
1949
1950
1951 template <int dim>
1952 Tensor<1, dim>
1953 FourierCosineFunction<dim>::gradient(const Point<dim> &p,
1954 const unsigned int component) const
1955 {
1956 (void)component;
1957 AssertIndexRange(component, 1);
1958 return -fourier_coefficients * std::sin(fourier_coefficients * p);
1959 }
1960
1961
1962
1963 template <int dim>
1964 double
1965 FourierCosineFunction<dim>::laplacian(const Point<dim> &p,
1966 const unsigned int component) const
1967 {
1968 (void)component;
1969 AssertIndexRange(component, 1);
1970 return (fourier_coefficients * fourier_coefficients) *
1971 (-std::cos(fourier_coefficients * p));
1972 }
1973
1974
1975
1976 /* ---------------------- FourierSineFunction ----------------------- */
1977
1978
1979
1980 template <int dim>
1981 FourierSineFunction<dim>::FourierSineFunction(
1982 const Tensor<1, dim> &fourier_coefficients)
1983 : Function<dim>(1)
1984 , fourier_coefficients(fourier_coefficients)
1985 {}
1986
1987
1988
1989 template <int dim>
1990 double
1991 FourierSineFunction<dim>::value(const Point<dim> &p,
1992 const unsigned int component) const
1993 {
1994 (void)component;
1995 AssertIndexRange(component, 1);
1996 return std::sin(fourier_coefficients * p);
1997 }
1998
1999
2000
2001 template <int dim>
2002 Tensor<1, dim>
2003 FourierSineFunction<dim>::gradient(const Point<dim> &p,
2004 const unsigned int component) const
2005 {
2006 (void)component;
2007 AssertIndexRange(component, 1);
2008 return fourier_coefficients * std::cos(fourier_coefficients * p);
2009 }
2010
2011
2012
2013 template <int dim>
2014 double
2015 FourierSineFunction<dim>::laplacian(const Point<dim> &p,
2016 const unsigned int component) const
2017 {
2018 (void)component;
2019 AssertIndexRange(component, 1);
2020 return (fourier_coefficients * fourier_coefficients) *
2021 (-std::sin(fourier_coefficients * p));
2022 }
2023
2024
2025
2026 /* ---------------------- FourierSineSum ----------------------- */
2027
2028
2029
2030 template <int dim>
2031 FourierSineSum<dim>::FourierSineSum(
2032 const std::vector<Point<dim>> &fourier_coefficients,
2033 const std::vector<double> &weights)
2034 : Function<dim>(1)
2035 , fourier_coefficients(fourier_coefficients)
2036 , weights(weights)
2037 {
2038 Assert(fourier_coefficients.size() > 0, ExcZero());
2039 Assert(fourier_coefficients.size() == weights.size(),
2040 ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2041 }
2042
2043
2044
2045 template <int dim>
2046 double
2047 FourierSineSum<dim>::value(const Point<dim> &p,
2048 const unsigned int component) const
2049 {
2050 (void)component;
2051 AssertIndexRange(component, 1);
2052
2053 const unsigned int n = weights.size();
2054 double sum = 0;
2055 for (unsigned int s = 0; s < n; ++s)
2056 sum += weights[s] * std::sin(fourier_coefficients[s] * p);
2057
2058 return sum;
2059 }
2060
2061
2062
2063 template <int dim>
2064 Tensor<1, dim>
2065 FourierSineSum<dim>::gradient(const Point<dim> &p,
2066 const unsigned int component) const
2067 {
2068 (void)component;
2069 AssertIndexRange(component, 1);
2070
2071 const unsigned int n = weights.size();
2072 Tensor<1, dim> sum;
2073 for (unsigned int s = 0; s < n; ++s)
2074 sum += fourier_coefficients[s] * std::cos(fourier_coefficients[s] * p);
2075
2076 return sum;
2077 }
2078
2079
2080
2081 template <int dim>
2082 double
2083 FourierSineSum<dim>::laplacian(const Point<dim> &p,
2084 const unsigned int component) const
2085 {
2086 (void)component;
2087 AssertIndexRange(component, 1);
2088
2089 const unsigned int n = weights.size();
2090 double sum = 0;
2091 for (unsigned int s = 0; s < n; ++s)
2092 sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2093 std::sin(fourier_coefficients[s] * p);
2094
2095 return sum;
2096 }
2097
2098
2099
2100 /* ---------------------- FourierCosineSum ----------------------- */
2101
2102
2103
2104 template <int dim>
2105 FourierCosineSum<dim>::FourierCosineSum(
2106 const std::vector<Point<dim>> &fourier_coefficients,
2107 const std::vector<double> &weights)
2108 : Function<dim>(1)
2109 , fourier_coefficients(fourier_coefficients)
2110 , weights(weights)
2111 {
2112 Assert(fourier_coefficients.size() > 0, ExcZero());
2113 Assert(fourier_coefficients.size() == weights.size(),
2114 ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2115 }
2116
2117
2118
2119 template <int dim>
2120 double
2121 FourierCosineSum<dim>::value(const Point<dim> &p,
2122 const unsigned int component) const
2123 {
2124 (void)component;
2125 AssertIndexRange(component, 1);
2126
2127 const unsigned int n = weights.size();
2128 double sum = 0;
2129 for (unsigned int s = 0; s < n; ++s)
2130 sum += weights[s] * std::cos(fourier_coefficients[s] * p);
2131
2132 return sum;
2133 }
2134
2135
2136
2137 template <int dim>
2138 Tensor<1, dim>
2139 FourierCosineSum<dim>::gradient(const Point<dim> &p,
2140 const unsigned int component) const
2141 {
2142 (void)component;
2143 AssertIndexRange(component, 1);
2144
2145 const unsigned int n = weights.size();
2146 Tensor<1, dim> sum;
2147 for (unsigned int s = 0; s < n; ++s)
2148 sum -= fourier_coefficients[s] * std::sin(fourier_coefficients[s] * p);
2149
2150 return sum;
2151 }
2152
2153
2154
2155 template <int dim>
2156 double
2157 FourierCosineSum<dim>::laplacian(const Point<dim> &p,
2158 const unsigned int component) const
2159 {
2160 (void)component;
2161 AssertIndexRange(component, 1);
2162
2163 const unsigned int n = weights.size();
2164 double sum = 0;
2165 for (unsigned int s = 0; s < n; ++s)
2166 sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2167 std::cos(fourier_coefficients[s] * p);
2168
2169 return sum;
2170 }
2171
2172
2173
2174 /* ---------------------- Monomial ----------------------- */
2175
2176
2177
2178 template <int dim, typename Number>
2179 Monomial<dim, Number>::Monomial(const Tensor<1, dim, Number> &exponents,
2180 const unsigned int n_components)
2181 : Function<dim, Number>(n_components)
2182 , exponents(exponents)
2183 {}
2184
2185
2186
2187 template <int dim, typename Number>
2188 Number
2189 Monomial<dim, Number>::value(const Point<dim> &p,
2190 const unsigned int component) const
2191 {
2192 (void)component;
2193 AssertIndexRange(component, this->n_components);
2194
2195 Number prod = 1;
2196 for (unsigned int s = 0; s < dim; ++s)
2197 {
2198 if (p[s] < 0)
2199 Assert(std::floor(exponents[s]) == exponents[s],
2200 ExcMessage("Exponentiation of a negative base number with "
2201 "a real exponent can't be performed."));
2202 prod *= std::pow(p[s], exponents[s]);
2203 }
2204 return prod;
2205 }
2206
2207
2208
2209 template <int dim, typename Number>
2210 void
2211 Monomial<dim, Number>::vector_value(const Point<dim> &p,
2212 Vector<Number> &values) const
2213 {
2214 Assert(values.size() == this->n_components,
2215 ExcDimensionMismatch(values.size(), this->n_components));
2216
2217 for (unsigned int i = 0; i < values.size(); ++i)
2218 values(i) = Monomial<dim, Number>::value(p, i);
2219 }
2220
2221
2222
2223 template <int dim, typename Number>
2224 Tensor<1, dim, Number>
2225 Monomial<dim, Number>::gradient(const Point<dim> &p,
2226 const unsigned int component) const
2227 {
2228 (void)component;
2229 AssertIndexRange(component, 1);
2230
2231 Tensor<1, dim> r;
2232 for (unsigned int d = 0; d < dim; ++d)
2233 {
2234 double prod = 1;
2235 for (unsigned int s = 0; s < dim; ++s)
2236 {
2237 if ((s == d) && (exponents[s] == 0) && (p[s] == 0))
2238 {
2239 prod = 0;
2240 break;
2241 }
2242 else
2243 {
2244 if (p[s] < 0)
2245 Assert(std::floor(exponents[s]) == exponents[s],
2246 ExcMessage(
2247 "Exponentiation of a negative base number with "
2248 "a real exponent can't be performed."));
2249 prod *=
2250 (s == d ? exponents[s] * std::pow(p[s], exponents[s] - 1) :
2251 std::pow(p[s], exponents[s]));
2252 }
2253 }
2254 r[d] = prod;
2255 }
2256
2257 return r;
2258 }
2259
2260
2261
2262 template <int dim, typename Number>
2263 void
2264 Monomial<dim, Number>::value_list(const std::vector<Point<dim>> &points,
2265 std::vector<Number> &values,
2266 const unsigned int component) const
2267 {
2268 Assert(values.size() == points.size(),
2269 ExcDimensionMismatch(values.size(), points.size()));
2270
2271 for (unsigned int i = 0; i < points.size(); ++i)
2272 values[i] = Monomial<dim, Number>::value(points[i], component);
2273 }
2274
2275
2276 template <int dim>
2277 Bessel1<dim>::Bessel1(const unsigned int order,
2278 const double wave_number,
2279 const Point<dim> center)
2280 : order(order)
2281 , wave_number(wave_number)
2282 , center(center)
2283 {
2284 Assert(wave_number >= 0., ExcMessage("wave_number must be nonnegative!"));
2285 }
2286
2287 template <int dim>
2288 double
2289 Bessel1<dim>::value(const Point<dim> &p, const unsigned int) const
2290 {
2291 Assert(dim == 2, ExcNotImplemented());
2292 const double r = p.distance(center);
2293 return std_cxx17::cyl_bessel_j(order, r * wave_number);
2294 }
2295
2296
2297 template <int dim>
2298 void
2299 Bessel1<dim>::value_list(const std::vector<Point<dim>> &points,
2300 std::vector<double> &values,
2301 const unsigned int) const
2302 {
2303 Assert(dim == 2, ExcNotImplemented());
2304 AssertDimension(points.size(), values.size());
2305 for (unsigned int k = 0; k < points.size(); ++k)
2306 {
2307 const double r = points[k].distance(center);
2308 values[k] = std_cxx17::cyl_bessel_j(order, r * wave_number);
2309 }
2310 }
2311
2312
2313 template <int dim>
2314 Tensor<1, dim>
2315 Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
2316 {
2317 Assert(dim == 2, ExcNotImplemented());
2318 const double r = p.distance(center);
2319 const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r;
2320 const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r;
2321
2322 const double dJn =
2323 (order == 0) ?
2324 (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2325 (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2326 std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2327 Tensor<1, dim> result;
2328 result[0] = wave_number * co * dJn;
2329 result[1] = wave_number * si * dJn;
2330 return result;
2331 }
2332
2333
2334
2335 template <int dim>
2336 void
2337 Bessel1<dim>::gradient_list(const std::vector<Point<dim>> &points,
2338 std::vector<Tensor<1, dim>> &gradients,
2339 const unsigned int) const
2340 {
2341 Assert(dim == 2, ExcNotImplemented());
2342 AssertDimension(points.size(), gradients.size());
2343 for (unsigned int k = 0; k < points.size(); ++k)
2344 {
2345 const Point<dim> &p = points[k];
2346 const double r = p.distance(center);
2347 const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r;
2348 const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r;
2349
2350 const double dJn =
2351 (order == 0) ?
2352 (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2353 (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2354 std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2355 Tensor<1, dim> &result = gradients[k];
2356 result[0] = wave_number * co * dJn;
2357 result[1] = wave_number * si * dJn;
2358 }
2359 }
2360
2361
2362
2363 namespace
2364 {
2365 // interpolate a data value from a table where ix denotes
2366 // the (lower) left endpoint of the interval to interpolate
2367 // in, and p_unit denotes the point in unit coordinates to do so.
2368 double
2369 interpolate(const Table<1, double> &data_values,
2370 const TableIndices<1> &ix,
2371 const Point<1> &xi)
2372 {
2373 return ((1 - xi[0]) * data_values[ix[0]] +
2374 xi[0] * data_values[ix[0] + 1]);
2375 }
2376
2377 double
2378 interpolate(const Table<2, double> &data_values,
2379 const TableIndices<2> &ix,
2380 const Point<2> &p_unit)
2381 {
2382 return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
2383 p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
2384 (1 - p_unit[1]) +
2385 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
2386 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
2387 p_unit[1]);
2388 }
2389
2390 double
2391 interpolate(const Table<3, double> &data_values,
2392 const TableIndices<3> &ix,
2393 const Point<3> &p_unit)
2394 {
2395 return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
2396 p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
2397 (1 - p_unit[1]) +
2398 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
2399 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
2400 p_unit[1]) *
2401 (1 - p_unit[2]) +
2402 (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
2403 p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
2404 (1 - p_unit[1]) +
2405 ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
2406 p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
2407 p_unit[1]) *
2408 p_unit[2]);
2409 }
2410
2411
2412 // Interpolate the gradient of a data value from a table where ix
2413 // denotes the lower left endpoint of the interval to interpolate
2414 // in, p_unit denotes the point in unit coordinates, and dx
2415 // denotes the width of the interval in each dimension.
2416 Tensor<1, 1>
2417 gradient_interpolate(const Table<1, double> &data_values,
2418 const TableIndices<1> &ix,
2419 const Point<1> &p_unit,
2420 const Point<1> &dx)
2421 {
2422 (void)p_unit;
2423 Tensor<1, 1> grad;
2424 grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
2425 return grad;
2426 }
2427
2428
2429 Tensor<1, 2>
2430 gradient_interpolate(const Table<2, double> &data_values,
2431 const TableIndices<2> &ix,
2432 const Point<2> &p_unit,
2433 const Point<2> &dx)
2434 {
2435 Tensor<1, 2> grad;
2436 double u00 = data_values[ix[0]][ix[1]],
2437 u01 = data_values[ix[0] + 1][ix[1]],
2438 u10 = data_values[ix[0]][ix[1] + 1],
2439 u11 = data_values[ix[0] + 1][ix[1] + 1];
2440
2441 grad[0] =
2442 ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
2443 grad[1] =
2444 ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
2445 return grad;
2446 }
2447
2448
2449 Tensor<1, 3>
2450 gradient_interpolate(const Table<3, double> &data_values,
2451 const TableIndices<3> &ix,
2452 const Point<3> &p_unit,
2453 const Point<3> &dx)
2454 {
2455 Tensor<1, 3> grad;
2456 double u000 = data_values[ix[0]][ix[1]][ix[2]],
2457 u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
2458 u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
2459 u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
2460 u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
2461 u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
2462 u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
2463 u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
2464
2465 grad[0] =
2466 ((1 - p_unit[2]) *
2467 ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
2468 p_unit[2] *
2469 ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
2470 dx[0];
2471 grad[1] =
2472 ((1 - p_unit[2]) *
2473 ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
2474 p_unit[2] *
2475 ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
2476 dx[1];
2477 grad[2] =
2478 ((1 - p_unit[1]) *
2479 ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
2480 p_unit[1] *
2481 ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
2482 dx[2];
2483
2484 return grad;
2485 }
2486 } // namespace
2487
2488
2489
2490 template <int dim>
2491 InterpolatedTensorProductGridData<dim>::InterpolatedTensorProductGridData(
2492 const std::array<std::vector<double>, dim> &coordinate_values,
2493 const Table<dim, double> &data_values)
2494 : coordinate_values(coordinate_values)
2495 , data_values(data_values)
2496 {
2497 for (unsigned int d = 0; d < dim; ++d)
2498 {
2499 Assert(
2500 coordinate_values[d].size() >= 2,
2501 ExcMessage(
2502 "Coordinate arrays must have at least two coordinate values!"));
2503 for (unsigned int i = 0; i < coordinate_values[d].size() - 1; ++i)
2504 Assert(
2505 coordinate_values[d][i] < coordinate_values[d][i + 1],
2506 ExcMessage(
2507 "Coordinate arrays must be sorted in strictly ascending order."));
2508
2509 Assert(data_values.size()[d] == coordinate_values[d].size(),
2510 ExcMessage(
2511 "Data and coordinate tables do not have the same size."));
2512 }
2513 }
2514
2515
2516
2517 template <int dim>
2518 InterpolatedTensorProductGridData<dim>::InterpolatedTensorProductGridData(
2519 std::array<std::vector<double>, dim> &&coordinate_values,
2520 Table<dim, double> &&data_values)
2521 : coordinate_values(std::move(coordinate_values))
2522 , data_values(std::move(data_values))
2523 {
2524 for (unsigned int d = 0; d < dim; ++d)
2525 {
2526 Assert(
2527 this->coordinate_values[d].size() >= 2,
2528 ExcMessage(
2529 "Coordinate arrays must have at least two coordinate values!"));
2530 for (unsigned int i = 0; i < this->coordinate_values[d].size() - 1; ++i)
2531 Assert(
2532 this->coordinate_values[d][i] < this->coordinate_values[d][i + 1],
2533 ExcMessage(
2534 "Coordinate arrays must be sorted in strictly ascending order."));
2535
2536 Assert(this->data_values.size()[d] == this->coordinate_values[d].size(),
2537 ExcMessage(
2538 "Data and coordinate tables do not have the same size."));
2539 }
2540 }
2541
2542
2543
2544 template <int dim>
2545 TableIndices<dim>
2546 InterpolatedTensorProductGridData<dim>::table_index_of_point(
2547 const Point<dim> &p) const
2548 {
2549 // find out where this data point lies, relative to the given
2550 // points. if we run all the way to the end of the range,
2551 // set the indices so that we will simply query the last of the
2552 // intervals, starting at x.size()-2 and going to x.size()-1.
2553 TableIndices<dim> ix;
2554 for (unsigned int d = 0; d < dim; ++d)
2555 {
2556 // get the index of the first element of the coordinate arrays that is
2557 // larger than p[d]
2558 ix[d] = (std::lower_bound(coordinate_values[d].begin(),
2559 coordinate_values[d].end(),
2560 p[d]) -
2561 coordinate_values[d].begin());
2562
2563 // the one we want is the index of the coordinate to the left, however,
2564 // so decrease it by one (unless we have a point to the left of all, in
2565 // which case we stay where we are; the formulas below are made in a way
2566 // that allow us to extend the function by a constant value)
2567 //
2568 // to make this work, if we got coordinate_values[d].end(), we actually
2569 // have to consider the last box which has index size()-2
2570 if (ix[d] == coordinate_values[d].size())
2571 ix[d] = coordinate_values[d].size() - 2;
2572 else if (ix[d] > 0)
2573 --ix[d];
2574 }
2575
2576 return ix;
2577 }
2578
2579
2580
2581 template <int dim>
2582 std::size_t
2583 InterpolatedTensorProductGridData<dim>::memory_consumption() const
2584 {
2585 return sizeof(*this) +
2586 MemoryConsumption::memory_consumption(coordinate_values) -
2587 sizeof(coordinate_values) +
2588 MemoryConsumption::memory_consumption(data_values) -
2589 sizeof(data_values);
2590 }
2591
2592
2593
2594 template <int dim>
2595 const Table<dim, double> &
2596 InterpolatedTensorProductGridData<dim>::get_data() const
2597 {
2598 return data_values;
2599 }
2600
2601
2602
2603 template <int dim>
2604 double
2605 InterpolatedTensorProductGridData<dim>::value(
2606 const Point<dim> &p,
2607 const unsigned int component) const
2608 {
2609 (void)component;
2610 Assert(
2611 component == 0,
2612 ExcMessage(
2613 "This is a scalar function object, the component can only be zero."));
2614
2615 // find the index in the data table of the cell containing the input point
2616 const TableIndices<dim> ix = table_index_of_point(p);
2617
2618 // now compute the relative point within the interval/rectangle/box
2619 // defined by the point coordinates found above. truncate below and
2620 // above to accommodate points that may lie outside the range
2621 Point<dim> p_unit;
2622 for (unsigned int d = 0; d < dim; ++d)
2623 p_unit[d] = std::max(std::min((p[d] - coordinate_values[d][ix[d]]) /
2624 (coordinate_values[d][ix[d] + 1] -
2625 coordinate_values[d][ix[d]]),
2626 1.),
2627 0.);
2628
2629 return interpolate(data_values, ix, p_unit);
2630 }
2631
2632
2633
2634 template <int dim>
2635 Tensor<1, dim>
2636 InterpolatedTensorProductGridData<dim>::gradient(
2637 const Point<dim> &p,
2638 const unsigned int component) const
2639 {
2640 (void)component;
2641 Assert(
2642 component == 0,
2643 ExcMessage(
2644 "This is a scalar function object, the component can only be zero."));
2645
2646 // find out where this data point lies
2647 const TableIndices<dim> ix = table_index_of_point(p);
2648
2649 Point<dim> dx;
2650 for (unsigned int d = 0; d < dim; ++d)
2651 dx[d] = coordinate_values[d][ix[d] + 1] - coordinate_values[d][ix[d]];
2652
2653 Point<dim> p_unit;
2654 for (unsigned int d = 0; d < dim; ++d)
2655 p_unit[d] =
2656 std::max(std::min((p[d] - coordinate_values[d][ix[d]]) / dx[d], 1.),
2657 0.0);
2658
2659 return gradient_interpolate(data_values, ix, p_unit, dx);
2660 }
2661
2662
2663
2664 template <int dim>
2665 InterpolatedUniformGridData<dim>::InterpolatedUniformGridData(
2666 const std::array<std::pair<double, double>, dim> &interval_endpoints,
2667 const std::array<unsigned int, dim> &n_subintervals,
2668 const Table<dim, double> &data_values)
2669 : interval_endpoints(interval_endpoints)
2670 , n_subintervals(n_subintervals)
2671 , data_values(data_values)
2672 {
2673 for (unsigned int d = 0; d < dim; ++d)
2674 {
2675 Assert(n_subintervals[d] >= 1,
2676 ExcMessage("There needs to be at least one subinterval in each "
2677 "coordinate direction."));
2678 Assert(interval_endpoints[d].first < interval_endpoints[d].second,
2679 ExcMessage("The interval in each coordinate direction needs "
2680 "to have positive size"));
2681 Assert(data_values.size()[d] == n_subintervals[d] + 1,
2682 ExcMessage("The data table does not have the correct size."));
2683 }
2684 }
2685
2686
2687
2688 template <int dim>
2689 InterpolatedUniformGridData<dim>::InterpolatedUniformGridData(
2690 std::array<std::pair<double, double>, dim> &&interval_endpoints,
2691 std::array<unsigned int, dim> &&n_subintervals,
2692 Table<dim, double> &&data_values)
2693 : interval_endpoints(std::move(interval_endpoints))
2694 , n_subintervals(std::move(n_subintervals))
2695 , data_values(std::move(data_values))
2696 {
2697 for (unsigned int d = 0; d < dim; ++d)
2698 {
2699 Assert(this->n_subintervals[d] >= 1,
2700 ExcMessage("There needs to be at least one subinterval in each "
2701 "coordinate direction."));
2702 Assert(this->interval_endpoints[d].first <
2703 this->interval_endpoints[d].second,
2704 ExcMessage("The interval in each coordinate direction needs "
2705 "to have positive size"));
2706 Assert(this->data_values.size()[d] == this->n_subintervals[d] + 1,
2707 ExcMessage("The data table does not have the correct size."));
2708 }
2709 }
2710
2711
2712
2713 template <int dim>
2714 double
2715 InterpolatedUniformGridData<dim>::value(const Point<dim> &p,
2716 const unsigned int component) const
2717 {
2718 (void)component;
2719 Assert(
2720 component == 0,
2721 ExcMessage(
2722 "This is a scalar function object, the component can only be zero."));
2723
2724 // find out where this data point lies, relative to the given
2725 // subdivision points
2726 TableIndices<dim> ix;
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 if (p[d] <= interval_endpoints[d].first)
2733 ix[d] = 0;
2734 else if (p[d] >= interval_endpoints[d].second - delta_x)
2735 ix[d] = n_subintervals[d] - 1;
2736 else
2737 ix[d] = static_cast<unsigned int>(
2738 (p[d] - interval_endpoints[d].first) / delta_x);
2739 }
2740
2741 // now compute the relative point within the interval/rectangle/box
2742 // defined by the point coordinates found above. truncate below and
2743 // above to accommodate points that may lie outside the range
2744 Point<dim> p_unit;
2745 for (unsigned int d = 0; d < dim; ++d)
2746 {
2747 const double delta_x =
2748 ((interval_endpoints[d].second - interval_endpoints[d].first) /
2749 n_subintervals[d]);
2750 p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
2751 ix[d] * delta_x) /
2752 delta_x,
2753 1.),
2754 0.);
2755 }
2756
2757 return interpolate(data_values, ix, p_unit);
2758 }
2759
2760
2761
2762 template <int dim>
2763 Tensor<1, dim>
2764 InterpolatedUniformGridData<dim>::gradient(const Point<dim> &p,
2765 const unsigned int component) const
2766 {
2767 (void)component;
2768 Assert(
2769 component == 0,
2770 ExcMessage(
2771 "This is a scalar function object, the component can only be zero."));
2772
2773 // find out where this data point lies, relative to the given
2774 // subdivision points
2775 TableIndices<dim> ix;
2776 for (unsigned int d = 0; d < dim; ++d)
2777 {
2778 const double delta_x = ((this->interval_endpoints[d].second -
2779 this->interval_endpoints[d].first) /
2780 this->n_subintervals[d]);
2781 if (p[d] <= this->interval_endpoints[d].first)
2782 ix[d] = 0;
2783 else if (p[d] >= this->interval_endpoints[d].second - delta_x)
2784 ix[d] = this->n_subintervals[d] - 1;
2785 else
2786 ix[d] = static_cast<unsigned int>(
2787 (p[d] - this->interval_endpoints[d].first) / delta_x);
2788 }
2789
2790 // now compute the relative point within the interval/rectangle/box
2791 // defined by the point coordinates found above. truncate below and
2792 // above to accommodate points that may lie outside the range
2793 Point<dim> p_unit;
2794 Point<dim> delta_x;
2795 for (unsigned int d = 0; d < dim; ++d)
2796 {
2797 delta_x[d] = ((this->interval_endpoints[d].second -
2798 this->interval_endpoints[d].first) /
2799 this->n_subintervals[d]);
2800 p_unit[d] =
2801 std::max(std::min((p[d] - this->interval_endpoints[d].first -
2802 ix[d] * delta_x[d]) /
2803 delta_x[d],
2804 1.),
2805 0.);
2806 }
2807
2808 return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
2809 }
2810
2811
2812
2813 template <int dim>
2814 std::size_t
2815 InterpolatedUniformGridData<dim>::memory_consumption() const
2816 {
2817 return sizeof(*this) + data_values.memory_consumption() -
2818 sizeof(data_values);
2819 }
2820
2821
2822
2823 template <int dim>
2824 const Table<dim, double> &
2825 InterpolatedUniformGridData<dim>::get_data() const
2826 {
2827 return data_values;
2828 }
2829
2830
2831
2832 /* ---------------------- Polynomial ----------------------- */
2833
2834
2835
2836 template <int dim>
2837 Polynomial<dim>::Polynomial(const Table<2, double> &exponents,
2838 const std::vector<double> &coefficients)
2839 : Function<dim>(1)
2840 , exponents(exponents)
2841 , coefficients(coefficients)
2842 {
2843 Assert(exponents.n_rows() == coefficients.size(),
2844 ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
2845 Assert(exponents.n_cols() == dim,
2846 ExcDimensionMismatch(exponents.n_cols(), dim));
2847 }
2848
2849
2850
2851 template <int dim>
2852 double
2853 Polynomial<dim>::value(const Point<dim> &p,
2854 const unsigned int component) const
2855 {
2856 (void)component;
2857 AssertIndexRange(component, 1);
2858
2859 double sum = 0;
2860 for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2861 {
2862 double prod = 1;
2863 for (unsigned int s = 0; s < dim; ++s)
2864 {
2865 if (p[s] < 0)
2866 Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
2867 ExcMessage("Exponentiation of a negative base number with "
2868 "a real exponent can't be performed."));
2869 prod *= std::pow(p[s], exponents[monom][s]);
2870 }
2871 sum += coefficients[monom] * prod;
2872 }
2873 return sum;
2874 }
2875
2876
2877
2878 template <int dim>
2879 void
2880 Polynomial<dim>::value_list(const std::vector<Point<dim>> &points,
2881 std::vector<double> &values,
2882 const unsigned int component) const
2883 {
2884 Assert(values.size() == points.size(),
2885 ExcDimensionMismatch(values.size(), points.size()));
2886
2887 for (unsigned int i = 0; i < points.size(); ++i)
2888 values[i] = Polynomial<dim>::value(points[i], component);
2889 }
2890
2891
2892
2893 template <int dim>
2894 Tensor<1, dim>
2895 Polynomial<dim>::gradient(const Point<dim> &p,
2896 const unsigned int component) const
2897 {
2898 (void)component;
2899 AssertIndexRange(component, 1);
2900
2901 Tensor<1, dim> r;
2902
2903 for (unsigned int d = 0; d < dim; ++d)
2904 {
2905 double sum = 0;
2906
2907 for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2908 {
2909 double prod = 1;
2910 for (unsigned int s = 0; s < dim; ++s)
2911 {
2912 if ((s == d) && (exponents[monom][s] == 0) && (p[s] == 0))
2913 {
2914 prod = 0;
2915 break;
2916 }
2917 else
2918 {
2919 if (p[s] < 0)
2920 Assert(std::floor(exponents[monom][s]) ==
2921 exponents[monom][s],
2922 ExcMessage(
2923 "Exponentiation of a negative base number with "
2924 "a real exponent can't be performed."));
2925 prod *=
2926 (s == d ? exponents[monom][s] *
2927 std::pow(p[s], exponents[monom][s] - 1) :
2928 std::pow(p[s], exponents[monom][s]));
2929 }
2930 }
2931 sum += coefficients[monom] * prod;
2932 }
2933 r[d] = sum;
2934 }
2935 return r;
2936 }
2937
2938
2939
2940 template <int dim>
2941 std::size_t
2942 Polynomial<dim>::memory_consumption() const
2943 {
2944 return sizeof(*this) + exponents.memory_consumption() - sizeof(exponents) +
2945 MemoryConsumption::memory_consumption(coefficients) -
2946 sizeof(coefficients);
2947 }
2948
2949 template <int dim>
2950 RayleighKotheVortex<dim>::RayleighKotheVortex(const double T)
2951 : Function<dim>(dim)
2952 , T(T)
2953 {
2954 AssertThrow(dim > 1, ExcNotImplemented());
2955 }
2956
2957
2958 template <int dim>
2959 void
2960 RayleighKotheVortex<dim>::vector_value(const Point<dim> &point,
2961 Vector<double> &values) const
2962 {
2963 const double pi_x = numbers::PI * point[0];
2964 const double pi_y = numbers::PI * point[1];
2965 const double pi_t = numbers::PI / T * this->get_time();
2966
2967 values[0] = -2 * std::cos(pi_t) *
2968 Utilities::fixed_power<2>(std::sin(pi_x)) * std::sin(pi_y) *
2969 std::cos(pi_y);
2970 values[1] = +2 * std::cos(pi_t) *
2971 Utilities::fixed_power<2>(std::sin(pi_y)) * std::sin(pi_x) *
2972 std::cos(pi_x);
2973
2974 if (dim == 3)
2975 values[2] = 0;
2976 }
2977
2978
2979 // explicit instantiations
2980 template class SquareFunction<1>;
2981 template class SquareFunction<2>;
2982 template class SquareFunction<3>;
2983 template class Q1WedgeFunction<1>;
2984 template class Q1WedgeFunction<2>;
2985 template class Q1WedgeFunction<3>;
2986 template class PillowFunction<1>;
2987 template class PillowFunction<2>;
2988 template class PillowFunction<3>;
2989 template class CosineFunction<1>;
2990 template class CosineFunction<2>;
2991 template class CosineFunction<3>;
2992 template class CosineGradFunction<1>;
2993 template class CosineGradFunction<2>;
2994 template class CosineGradFunction<3>;
2995 template class ExpFunction<1>;
2996 template class ExpFunction<2>;
2997 template class ExpFunction<3>;
2998 template class JumpFunction<1>;
2999 template class JumpFunction<2>;
3000 template class JumpFunction<3>;
3001 template class FourierCosineFunction<1>;
3002 template class FourierCosineFunction<2>;
3003 template class FourierCosineFunction<3>;
3004 template class FourierSineFunction<1>;
3005 template class FourierSineFunction<2>;
3006 template class FourierSineFunction<3>;
3007 template class FourierCosineSum<1>;
3008 template class FourierCosineSum<2>;
3009 template class FourierCosineSum<3>;
3010 template class FourierSineSum<1>;
3011 template class FourierSineSum<2>;
3012 template class FourierSineSum<3>;
3013 template class SlitSingularityFunction<2>;
3014 template class SlitSingularityFunction<3>;
3015 template class Monomial<1>;
3016 template class Monomial<2>;
3017 template class Monomial<3>;
3018 template class Monomial<1, float>;
3019 template class Monomial<2, float>;
3020 template class Monomial<3, float>;
3021 template class Bessel1<1>;
3022 template class Bessel1<2>;
3023 template class Bessel1<3>;
3024 template class InterpolatedTensorProductGridData<1>;
3025 template class InterpolatedTensorProductGridData<2>;
3026 template class InterpolatedTensorProductGridData<3>;
3027 template class InterpolatedUniformGridData<1>;
3028 template class InterpolatedUniformGridData<2>;
3029 template class InterpolatedUniformGridData<3>;
3030 template class Polynomial<1>;
3031 template class Polynomial<2>;
3032 template class Polynomial<3>;
3033 template class RayleighKotheVortex<1>;
3034 template class RayleighKotheVortex<2>;
3035 template class RayleighKotheVortex<3>;
3036} // namespace Functions
3037
3038DEAL_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:2315
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:2299
Functions::Bessel1::value
virtual double value(const Point< dim > &points, const unsigned int component=0) const override
Definition function_lib.cc:2289
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:2337
Functions::Bessel1::Bessel1
Bessel1(const unsigned int order, const double wave_number, const Point< dim > center=Point< dim >())
Definition function_lib.cc:2277
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:1965
Functions::FourierCosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1953
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:2157
Functions::FourierCosineSum::FourierCosineSum
FourierCosineSum(const std::vector< Point< dim > > &fourier_coefficients, const std::vector< double > &weights)
Definition function_lib.cc:2105
Functions::FourierCosineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2139
Functions::FourierCosineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2121
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:1981
Functions::FourierSineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2015
Functions::FourierSineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:1991
Functions::FourierSineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2003
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:2065
Functions::FourierSineSum::FourierSineSum
FourierSineSum(const std::vector< Point< dim > > &fourier_coefficients, const std::vector< double > &weights)
Definition function_lib.cc:2031
Functions::FourierSineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2047
Functions::FourierSineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2083
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:2605
Functions::InterpolatedTensorProductGridData::InterpolatedTensorProductGridData
InterpolatedTensorProductGridData(const std::array< std::vector< double >, dim > &coordinate_values, const Table< dim, double > &data_values)
Definition function_lib.cc:2491
Functions::InterpolatedTensorProductGridData::get_data
const Table< dim, double > & get_data() const
Definition function_lib.cc:2596
Functions::InterpolatedTensorProductGridData::table_index_of_point
TableIndices< dim > table_index_of_point(const Point< dim > &p) const
Definition function_lib.cc:2546
Functions::InterpolatedTensorProductGridData::memory_consumption
virtual std::size_t memory_consumption() const override
Definition function_lib.cc:2583
Functions::InterpolatedTensorProductGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2636
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:2825
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:2764
Functions::InterpolatedUniformGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2715
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:2665
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:2815
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:2189
Functions::Monomial::gradient
virtual Tensor< 1, dim, Number > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2225
Functions::Monomial::Monomial
Monomial(const Tensor< 1, dim, Number > &exponents, const unsigned int n_components=1)
Definition function_lib.cc:2179
Functions::Monomial::vector_value
virtual void vector_value(const Point< dim > &p, Vector< Number > &values) const override
Definition function_lib.cc:2211
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:2264
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:2880
Functions::Polynomial::memory_consumption
virtual std::size_t memory_consumption() const override
Definition function_lib.cc:2942
Functions::Polynomial::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2853
Functions::Polynomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition function_lib.cc:2895
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:2837
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:2950
Functions::RayleighKotheVortex::vector_value
virtual void vector_value(const Point< dim > &point, Vector< double > &values) const override
Definition function_lib.cc:2960
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:503
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
center
Point< 3 > center
Definition data_out_base.cc:267
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:1638
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition exceptions.h:2008
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2053
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:1739
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1814
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
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
numbers::PI_2
static constexpr double PI_2
Definition numbers.h:264
numbers::PI
static constexpr double PI
Definition numbers.h:259
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.11.0