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