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