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