Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
function_lib.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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 
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 
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 
423 
424  template <int dim>
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 
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 
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 
1170 
1171 
1172  double
1173  LSingularityFunction::value(const Point<2> &p, const unsigned int) const
1174  {
1175  double x = p(0);
1176  double y = p(1);
1177 
1178  if ((x >= 0) && (y >= 0))
1179  return 0.;
1180 
1181  double phi = std::atan2(y, -x) + numbers::PI;
1182  double r2 = x * x + y * y;
1183 
1184  return std::pow(r2, 1. / 3.) * std::sin(2. / 3. * phi);
1185  }
1186 
1187 
1188  void
1189  LSingularityFunction::value_list(const std::vector<Point<2>> &points,
1190  std::vector<double> & values,
1191  const unsigned int) const
1192  {
1193  Assert(values.size() == points.size(),
1194  ExcDimensionMismatch(values.size(), points.size()));
1195 
1196  for (unsigned int i = 0; i < points.size(); ++i)
1197  {
1198  double x = points[i](0);
1199  double y = points[i](1);
1200 
1201  if ((x >= 0) && (y >= 0))
1202  values[i] = 0.;
1203  else
1204  {
1205  double phi = std::atan2(y, -x) + numbers::PI;
1206  double r2 = x * x + y * y;
1207 
1208  values[i] = std::pow(r2, 1. / 3.) * std::sin(2. / 3. * phi);
1209  }
1210  }
1211  }
1212 
1213 
1214  void
1216  const std::vector<Point<2>> &points,
1217  std::vector<Vector<double>> &values) const
1218  {
1219  Assert(values.size() == points.size(),
1220  ExcDimensionMismatch(values.size(), points.size()));
1221 
1222  for (unsigned int i = 0; i < points.size(); ++i)
1223  {
1224  Assert(values[i].size() == 1,
1225  ExcDimensionMismatch(values[i].size(), 1));
1226  double x = points[i](0);
1227  double y = points[i](1);
1228 
1229  if ((x >= 0) && (y >= 0))
1230  values[i](0) = 0.;
1231  else
1232  {
1233  double phi = std::atan2(y, -x) + numbers::PI;
1234  double r2 = x * x + y * y;
1235 
1236  values[i](0) = std::pow(r2, 1. / 3.) * std::sin(2. / 3. * phi);
1237  }
1238  }
1239  }
1240 
1241 
1242  double
1243  LSingularityFunction::laplacian(const Point<2> &, const unsigned int) const
1244  {
1245  return 0.;
1246  }
1247 
1248 
1249  void
1251  std::vector<double> & values,
1252  const unsigned int) const
1253  {
1254  Assert(values.size() == points.size(),
1255  ExcDimensionMismatch(values.size(), points.size()));
1256 
1257  for (unsigned int i = 0; i < points.size(); ++i)
1258  values[i] = 0.;
1259  }
1260 
1261 
1262  Tensor<1, 2>
1263  LSingularityFunction::gradient(const Point<2> &p, const unsigned int) const
1264  {
1265  double x = p(0);
1266  double y = p(1);
1267  double phi = std::atan2(y, -x) + numbers::PI;
1268  double r43 = std::pow(x * x + y * y, 2. / 3.);
1269 
1270  Tensor<1, 2> result;
1271  result[0] = 2. / 3. *
1272  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) /
1273  r43;
1274  result[1] = 2. / 3. *
1275  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) /
1276  r43;
1277  return result;
1278  }
1279 
1280 
1281  void
1282  LSingularityFunction::gradient_list(const std::vector<Point<2>> &points,
1283  std::vector<Tensor<1, 2>> & gradients,
1284  const unsigned int) const
1285  {
1286  Assert(gradients.size() == points.size(),
1287  ExcDimensionMismatch(gradients.size(), points.size()));
1288 
1289  for (unsigned int i = 0; i < points.size(); ++i)
1290  {
1291  const Point<2> &p = points[i];
1292  double x = p(0);
1293  double y = p(1);
1294  double phi = std::atan2(y, -x) + numbers::PI;
1295  double r43 = std::pow(x * x + y * y, 2. / 3.);
1296 
1297  gradients[i][0] =
1298  2. / 3. *
1299  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1300  gradients[i][1] =
1301  2. / 3. *
1302  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1303  }
1304  }
1305 
1306 
1307  void
1309  const std::vector<Point<2>> & points,
1310  std::vector<std::vector<Tensor<1, 2>>> &gradients) const
1311  {
1312  Assert(gradients.size() == points.size(),
1313  ExcDimensionMismatch(gradients.size(), points.size()));
1314 
1315  for (unsigned int i = 0; i < points.size(); ++i)
1316  {
1317  Assert(gradients[i].size() == 1,
1318  ExcDimensionMismatch(gradients[i].size(), 1));
1319  const Point<2> &p = points[i];
1320  double x = p(0);
1321  double y = p(1);
1322  double phi = std::atan2(y, -x) + numbers::PI;
1323  double r43 = std::pow(x * x + y * y, 2. / 3.);
1324 
1325  gradients[i][0][0] =
1326  2. / 3. *
1327  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1328  gradients[i][0][1] =
1329  2. / 3. *
1330  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1331  }
1332  }
1333 
1335 
1337  : Function<2>(2)
1338  {}
1339 
1340 
1341  double
1342  LSingularityGradFunction::value(const Point<2> &p, const unsigned int d) const
1343  {
1344  AssertIndexRange(d, 2);
1345 
1346  const double x = p(0);
1347  const double y = p(1);
1348  const double phi = std::atan2(y, -x) + numbers::PI;
1349  const double r43 = std::pow(x * x + y * y, 2. / 3.);
1350 
1351  return 2. / 3. *
1352  (std::sin(2. / 3. * phi) * p(d) +
1353  (d == 0 ? (std::cos(2. / 3. * phi) * y) :
1354  (-std::cos(2. / 3. * phi) * x))) /
1355  r43;
1356  }
1357 
1358 
1359  void
1361  std::vector<double> & values,
1362  const unsigned int d) const
1363  {
1364  AssertIndexRange(d, 2);
1365  AssertDimension(values.size(), points.size());
1366 
1367  for (unsigned int i = 0; i < points.size(); ++i)
1368  {
1369  const Point<2> &p = points[i];
1370  const double x = p(0);
1371  const double y = p(1);
1372  const double phi = std::atan2(y, -x) + numbers::PI;
1373  const double r43 = std::pow(x * x + y * y, 2. / 3.);
1374 
1375  values[i] = 2. / 3. *
1376  (std::sin(2. / 3. * phi) * p(d) +
1377  (d == 0 ? (std::cos(2. / 3. * phi) * y) :
1378  (-std::cos(2. / 3. * phi) * x))) /
1379  r43;
1380  }
1381  }
1382 
1383 
1384  void
1386  const std::vector<Point<2>> &points,
1387  std::vector<Vector<double>> &values) const
1388  {
1389  Assert(values.size() == points.size(),
1390  ExcDimensionMismatch(values.size(), points.size()));
1391 
1392  for (unsigned int i = 0; i < points.size(); ++i)
1393  {
1394  AssertDimension(values[i].size(), 2);
1395  const Point<2> &p = points[i];
1396  const double x = p(0);
1397  const double y = p(1);
1398  const double phi = std::atan2(y, -x) + numbers::PI;
1399  const double r43 = std::pow(x * x + y * y, 2. / 3.);
1400 
1401  values[i](0) =
1402  2. / 3. *
1403  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1404  values[i](1) =
1405  2. / 3. *
1406  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1407  }
1408  }
1409 
1410 
1411  double
1413  const unsigned int) const
1414  {
1415  return 0.;
1416  }
1417 
1418 
1419  void
1421  std::vector<double> & values,
1422  const unsigned int) const
1423  {
1424  Assert(values.size() == points.size(),
1425  ExcDimensionMismatch(values.size(), points.size()));
1426 
1427  for (unsigned int i = 0; i < points.size(); ++i)
1428  values[i] = 0.;
1429  }
1430 
1431 
1432 
1433  Tensor<1, 2>
1435  const unsigned int /*component*/) const
1436  {
1437  Assert(false, ExcNotImplemented());
1438  return {};
1439  }
1440 
1441 
1442  void
1444  const std::vector<Point<2>> & /*points*/,
1445  std::vector<Tensor<1, 2>> & /*gradients*/,
1446  const unsigned int /*component*/) const
1447  {
1448  Assert(false, ExcNotImplemented());
1449  }
1450 
1451 
1452  void
1454  const std::vector<Point<2>> & /*points*/,
1455  std::vector<std::vector<Tensor<1, 2>>> & /*gradients*/) const
1456  {
1457  Assert(false, ExcNotImplemented());
1458  }
1459 
1461 
1462  template <int dim>
1463  double
1465  const unsigned int) const
1466  {
1467  double x = p(0);
1468  double y = p(1);
1469 
1470  double phi = std::atan2(x, y) + numbers::PI;
1471  double r2 = x * x + y * y;
1472 
1473  return std::pow(r2, .25) * std::sin(.5 * phi);
1474  }
1475 
1476 
1477  template <int dim>
1478  void
1480  const std::vector<Point<dim>> &points,
1481  std::vector<double> & values,
1482  const unsigned int) const
1483  {
1484  Assert(values.size() == points.size(),
1485  ExcDimensionMismatch(values.size(), points.size()));
1486 
1487  for (unsigned int i = 0; i < points.size(); ++i)
1488  {
1489  double x = points[i](0);
1490  double y = points[i](1);
1491 
1492  double phi = std::atan2(x, y) + numbers::PI;
1493  double r2 = x * x + y * y;
1494 
1495  values[i] = std::pow(r2, .25) * std::sin(.5 * phi);
1496  }
1497  }
1498 
1499 
1500  template <int dim>
1501  void
1503  const std::vector<Point<dim>> &points,
1504  std::vector<Vector<double>> & values) const
1505  {
1506  Assert(values.size() == points.size(),
1507  ExcDimensionMismatch(values.size(), points.size()));
1508 
1509  for (unsigned int i = 0; i < points.size(); ++i)
1510  {
1511  Assert(values[i].size() == 1,
1512  ExcDimensionMismatch(values[i].size(), 1));
1513 
1514  double x = points[i](0);
1515  double y = points[i](1);
1516 
1517  double phi = std::atan2(x, y) + numbers::PI;
1518  double r2 = x * x + y * y;
1519 
1520  values[i](0) = std::pow(r2, .25) * std::sin(.5 * phi);
1521  }
1522  }
1523 
1524 
1525  template <int dim>
1526  double
1528  const unsigned int) const
1529  {
1530  return 0.;
1531  }
1532 
1533 
1534  template <int dim>
1535  void
1537  const std::vector<Point<dim>> &points,
1538  std::vector<double> & values,
1539  const unsigned int) const
1540  {
1541  Assert(values.size() == points.size(),
1542  ExcDimensionMismatch(values.size(), points.size()));
1543 
1544  for (unsigned int i = 0; i < points.size(); ++i)
1545  values[i] = 0.;
1546  }
1547 
1548 
1549  template <int dim>
1552  const unsigned int) const
1553  {
1554  double x = p(0);
1555  double y = p(1);
1556  double phi = std::atan2(x, y) + numbers::PI;
1557  double r64 = std::pow(x * x + y * y, 3. / 4.);
1558 
1559  Tensor<1, dim> result;
1560  result[0] = 1. / 2. *
1561  (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) /
1562  r64;
1563  result[1] = 1. / 2. *
1564  (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) /
1565  r64;
1566  return result;
1567  }
1568 
1569 
1570  template <int dim>
1571  void
1573  const std::vector<Point<dim>> &points,
1574  std::vector<Tensor<1, dim>> & gradients,
1575  const unsigned int) const
1576  {
1577  Assert(gradients.size() == points.size(),
1578  ExcDimensionMismatch(gradients.size(), points.size()));
1579 
1580  for (unsigned int i = 0; i < points.size(); ++i)
1581  {
1582  const Point<dim> &p = points[i];
1583  double x = p(0);
1584  double y = p(1);
1585  double phi = std::atan2(x, y) + numbers::PI;
1586  double r64 = std::pow(x * x + y * y, 3. / 4.);
1587 
1588  gradients[i][0] =
1589  1. / 2. *
1590  (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) / r64;
1591  gradients[i][1] =
1592  1. / 2. *
1593  (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) / r64;
1594  for (unsigned int d = 2; d < dim; ++d)
1595  gradients[i][d] = 0.;
1596  }
1597  }
1598 
1599  template <int dim>
1600  void
1602  const std::vector<Point<dim>> & points,
1603  std::vector<std::vector<Tensor<1, dim>>> &gradients) const
1604  {
1605  Assert(gradients.size() == points.size(),
1606  ExcDimensionMismatch(gradients.size(), points.size()));
1607 
1608  for (unsigned int i = 0; i < points.size(); ++i)
1609  {
1610  Assert(gradients[i].size() == 1,
1611  ExcDimensionMismatch(gradients[i].size(), 1));
1612 
1613  const Point<dim> &p = points[i];
1614  double x = p(0);
1615  double y = p(1);
1616  double phi = std::atan2(x, y) + numbers::PI;
1617  double r64 = std::pow(x * x + y * y, 3. / 4.);
1618 
1619  gradients[i][0][0] =
1620  1. / 2. *
1621  (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) / r64;
1622  gradients[i][0][1] =
1623  1. / 2. *
1624  (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) / r64;
1625  for (unsigned int d = 2; d < dim; ++d)
1626  gradients[i][0][d] = 0.;
1627  }
1628  }
1629 
1631 
1632 
1633  double
1635  const unsigned int) const
1636  {
1637  double x = p(0);
1638  double y = p(1);
1639 
1640  double phi = std::atan2(x, y) + numbers::PI;
1641  double r2 = x * x + y * y;
1642 
1643  return std::pow(r2, .125) * std::sin(.25 * phi);
1644  }
1645 
1646 
1647  void
1649  std::vector<double> & values,
1650  const unsigned int) const
1651  {
1652  Assert(values.size() == points.size(),
1653  ExcDimensionMismatch(values.size(), points.size()));
1654 
1655  for (unsigned int i = 0; i < points.size(); ++i)
1656  {
1657  double x = points[i](0);
1658  double y = points[i](1);
1659 
1660  double phi = std::atan2(x, y) + numbers::PI;
1661  double r2 = x * x + y * y;
1662 
1663  values[i] = std::pow(r2, .125) * std::sin(.25 * phi);
1664  }
1665  }
1666 
1667 
1668  void
1670  const std::vector<Point<2>> &points,
1671  std::vector<Vector<double>> &values) const
1672  {
1673  Assert(values.size() == points.size(),
1674  ExcDimensionMismatch(values.size(), points.size()));
1675 
1676  for (unsigned int i = 0; i < points.size(); ++i)
1677  {
1678  Assert(values[i].size() == 1,
1679  ExcDimensionMismatch(values[i].size(), 1));
1680 
1681  double x = points[i](0);
1682  double y = points[i](1);
1683 
1684  double phi = std::atan2(x, y) + numbers::PI;
1685  double r2 = x * x + y * y;
1686 
1687  values[i](0) = std::pow(r2, .125) * std::sin(.25 * phi);
1688  }
1689  }
1690 
1691 
1692  double
1694  const unsigned int) const
1695  {
1696  return 0.;
1697  }
1698 
1699 
1700  void
1702  const std::vector<Point<2>> &points,
1703  std::vector<double> & values,
1704  const unsigned int) const
1705  {
1706  Assert(values.size() == points.size(),
1707  ExcDimensionMismatch(values.size(), points.size()));
1708 
1709  for (unsigned int i = 0; i < points.size(); ++i)
1710  values[i] = 0.;
1711  }
1712 
1713 
1714  Tensor<1, 2>
1716  const unsigned int) const
1717  {
1718  double x = p(0);
1719  double y = p(1);
1720  double phi = std::atan2(x, y) + numbers::PI;
1721  double r78 = std::pow(x * x + y * y, 7. / 8.);
1722 
1723 
1724  Tensor<1, 2> result;
1725  result[0] = 1. / 4. *
1726  (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) /
1727  r78;
1728  result[1] = 1. / 4. *
1729  (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) /
1730  r78;
1731  return result;
1732  }
1733 
1734 
1735  void
1737  const std::vector<Point<2>> &points,
1738  std::vector<Tensor<1, 2>> & gradients,
1739  const unsigned int) const
1740  {
1741  Assert(gradients.size() == points.size(),
1742  ExcDimensionMismatch(gradients.size(), points.size()));
1743 
1744  for (unsigned int i = 0; i < points.size(); ++i)
1745  {
1746  const Point<2> &p = points[i];
1747  double x = p(0);
1748  double y = p(1);
1749  double phi = std::atan2(x, y) + numbers::PI;
1750  double r78 = std::pow(x * x + y * y, 7. / 8.);
1751 
1752  gradients[i][0] =
1753  1. / 4. *
1754  (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) / r78;
1755  gradients[i][1] =
1756  1. / 4. *
1757  (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) / r78;
1758  }
1759  }
1760 
1761 
1762  void
1764  const std::vector<Point<2>> & points,
1765  std::vector<std::vector<Tensor<1, 2>>> &gradients) const
1766  {
1767  Assert(gradients.size() == points.size(),
1768  ExcDimensionMismatch(gradients.size(), points.size()));
1769 
1770  for (unsigned int i = 0; i < points.size(); ++i)
1771  {
1772  Assert(gradients[i].size() == 1,
1773  ExcDimensionMismatch(gradients[i].size(), 1));
1774 
1775  const Point<2> &p = points[i];
1776  double x = p(0);
1777  double y = p(1);
1778  double phi = std::atan2(x, y) + numbers::PI;
1779  double r78 = std::pow(x * x + y * y, 7. / 8.);
1780 
1781  gradients[i][0][0] =
1782  1. / 4. *
1783  (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) / r78;
1784  gradients[i][0][1] =
1785  1. / 4. *
1786  (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) / r78;
1787  }
1788  }
1789 
1791 
1792  template <int dim>
1794  const double steepness)
1795  : direction(direction)
1796  , steepness(steepness)
1797  {
1798  switch (dim)
1799  {
1800  case 1:
1801  angle = 0;
1802  break;
1803  case 2:
1804  angle = std::atan2(direction(0), direction(1));
1805  break;
1806  case 3:
1807  Assert(false, ExcNotImplemented());
1808  }
1809  sine = std::sin(angle);
1810  cosine = std::cos(angle);
1811  }
1812 
1813 
1814 
1815  template <int dim>
1816  double
1817  JumpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1818  {
1819  double x = steepness * (-cosine * p(0) + sine * p(1));
1820  return -std::atan(x);
1821  }
1822 
1823 
1824 
1825  template <int dim>
1826  void
1828  std::vector<double> & values,
1829  const unsigned int) const
1830  {
1831  Assert(values.size() == p.size(),
1832  ExcDimensionMismatch(values.size(), p.size()));
1833 
1834  for (unsigned int i = 0; i < p.size(); ++i)
1835  {
1836  double x = steepness * (-cosine * p[i](0) + sine * p[i](1));
1837  values[i] = -std::atan(x);
1838  }
1839  }
1840 
1841 
1842  template <int dim>
1843  double
1844  JumpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1845  {
1846  double x = steepness * (-cosine * p(0) + sine * p(1));
1847  double r = 1 + x * x;
1848  return 2 * steepness * steepness * x / (r * r);
1849  }
1850 
1851 
1852  template <int dim>
1853  void
1855  std::vector<double> & values,
1856  const unsigned int) const
1857  {
1858  Assert(values.size() == p.size(),
1859  ExcDimensionMismatch(values.size(), p.size()));
1860 
1861  double f = 2 * steepness * steepness;
1862 
1863  for (unsigned int i = 0; i < p.size(); ++i)
1864  {
1865  double x = steepness * (-cosine * p[i](0) + sine * p[i](1));
1866  double r = 1 + x * x;
1867  values[i] = f * x / (r * r);
1868  }
1869  }
1870 
1871 
1872 
1873  template <int dim>
1875  JumpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1876  {
1877  double x = steepness * (-cosine * p(0) + sine * p(1));
1878  double r = -steepness * (1 + x * x);
1879  Tensor<1, dim> erg;
1880  erg[0] = cosine * r;
1881  erg[1] = sine * r;
1882  return erg;
1883  }
1884 
1885 
1886 
1887  template <int dim>
1888  void
1890  std::vector<Tensor<1, dim>> & gradients,
1891  const unsigned int) const
1892  {
1893  Assert(gradients.size() == p.size(),
1894  ExcDimensionMismatch(gradients.size(), p.size()));
1895 
1896  for (unsigned int i = 0; i < p.size(); ++i)
1897  {
1898  double x = steepness * (cosine * p[i](0) + sine * p[i](1));
1899  double r = -steepness * (1 + x * x);
1900  gradients[i][0] = cosine * r;
1901  gradients[i][1] = sine * r;
1902  }
1903  }
1904 
1905 
1906 
1907  template <int dim>
1908  std::size_t
1910  {
1911  // only simple data elements, so
1912  // use sizeof operator
1913  return sizeof(*this);
1914  }
1915 
1916 
1917 
1918  /* ---------------------- FourierCosineFunction ----------------------- */
1919 
1920 
1921  template <int dim>
1923  const Tensor<1, dim> &fourier_coefficients)
1924  : Function<dim>(1)
1925  , fourier_coefficients(fourier_coefficients)
1926  {}
1927 
1928 
1929 
1930  template <int dim>
1931  double
1933  const unsigned int component) const
1934  {
1935  (void)component;
1936  AssertIndexRange(component, 1);
1937  return std::cos(fourier_coefficients * p);
1938  }
1939 
1940 
1941 
1942  template <int dim>
1945  const unsigned int component) const
1946  {
1947  (void)component;
1948  AssertIndexRange(component, 1);
1949  return -fourier_coefficients * std::sin(fourier_coefficients * p);
1950  }
1951 
1952 
1953 
1954  template <int dim>
1955  double
1957  const unsigned int component) const
1958  {
1959  (void)component;
1960  AssertIndexRange(component, 1);
1961  return (fourier_coefficients * fourier_coefficients) *
1962  (-std::cos(fourier_coefficients * p));
1963  }
1964 
1965 
1966 
1967  /* ---------------------- FourierSineFunction ----------------------- */
1968 
1969 
1970 
1971  template <int dim>
1973  const Tensor<1, dim> &fourier_coefficients)
1974  : Function<dim>(1)
1975  , fourier_coefficients(fourier_coefficients)
1976  {}
1977 
1978 
1979 
1980  template <int dim>
1981  double
1983  const unsigned int component) const
1984  {
1985  (void)component;
1986  AssertIndexRange(component, 1);
1987  return std::sin(fourier_coefficients * p);
1988  }
1989 
1990 
1991 
1992  template <int dim>
1995  const unsigned int component) const
1996  {
1997  (void)component;
1998  AssertIndexRange(component, 1);
1999  return fourier_coefficients * std::cos(fourier_coefficients * p);
2000  }
2001 
2002 
2003 
2004  template <int dim>
2005  double
2007  const unsigned int component) const
2008  {
2009  (void)component;
2010  AssertIndexRange(component, 1);
2011  return (fourier_coefficients * fourier_coefficients) *
2012  (-std::sin(fourier_coefficients * p));
2013  }
2014 
2015 
2016 
2017  /* ---------------------- FourierSineSum ----------------------- */
2018 
2019 
2020 
2021  template <int dim>
2023  const std::vector<Point<dim>> &fourier_coefficients,
2024  const std::vector<double> & weights)
2025  : Function<dim>(1)
2026  , fourier_coefficients(fourier_coefficients)
2027  , weights(weights)
2028  {
2029  Assert(fourier_coefficients.size() > 0, ExcZero());
2030  Assert(fourier_coefficients.size() == weights.size(),
2032  }
2033 
2034 
2035 
2036  template <int dim>
2037  double
2039  const unsigned int component) const
2040  {
2041  (void)component;
2042  AssertIndexRange(component, 1);
2043 
2044  const unsigned int n = weights.size();
2045  double sum = 0;
2046  for (unsigned int s = 0; s < n; ++s)
2047  sum += weights[s] * std::sin(fourier_coefficients[s] * p);
2048 
2049  return sum;
2050  }
2051 
2052 
2053 
2054  template <int dim>
2057  const unsigned int component) const
2058  {
2059  (void)component;
2060  AssertIndexRange(component, 1);
2061 
2062  const unsigned int n = weights.size();
2064  for (unsigned int s = 0; s < n; ++s)
2065  sum += fourier_coefficients[s] * std::cos(fourier_coefficients[s] * p);
2066 
2067  return sum;
2068  }
2069 
2070 
2071 
2072  template <int dim>
2073  double
2075  const unsigned int component) const
2076  {
2077  (void)component;
2078  AssertIndexRange(component, 1);
2079 
2080  const unsigned int n = weights.size();
2081  double sum = 0;
2082  for (unsigned int s = 0; s < n; ++s)
2083  sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2084  std::sin(fourier_coefficients[s] * p);
2085 
2086  return sum;
2087  }
2088 
2089 
2090 
2091  /* ---------------------- FourierCosineSum ----------------------- */
2092 
2093 
2094 
2095  template <int dim>
2097  const std::vector<Point<dim>> &fourier_coefficients,
2098  const std::vector<double> & weights)
2099  : Function<dim>(1)
2100  , fourier_coefficients(fourier_coefficients)
2101  , weights(weights)
2102  {
2103  Assert(fourier_coefficients.size() > 0, ExcZero());
2104  Assert(fourier_coefficients.size() == weights.size(),
2106  }
2107 
2108 
2109 
2110  template <int dim>
2111  double
2113  const unsigned int component) const
2114  {
2115  (void)component;
2116  AssertIndexRange(component, 1);
2117 
2118  const unsigned int n = weights.size();
2119  double sum = 0;
2120  for (unsigned int s = 0; s < n; ++s)
2121  sum += weights[s] * std::cos(fourier_coefficients[s] * p);
2122 
2123  return sum;
2124  }
2125 
2126 
2127 
2128  template <int dim>
2131  const unsigned int component) const
2132  {
2133  (void)component;
2134  AssertIndexRange(component, 1);
2135 
2136  const unsigned int n = weights.size();
2138  for (unsigned int s = 0; s < n; ++s)
2139  sum -= fourier_coefficients[s] * std::sin(fourier_coefficients[s] * p);
2140 
2141  return sum;
2142  }
2143 
2144 
2145 
2146  template <int dim>
2147  double
2149  const unsigned int component) const
2150  {
2151  (void)component;
2152  AssertIndexRange(component, 1);
2153 
2154  const unsigned int n = weights.size();
2155  double sum = 0;
2156  for (unsigned int s = 0; s < n; ++s)
2157  sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2158  std::cos(fourier_coefficients[s] * p);
2159 
2160  return sum;
2161  }
2162 
2163 
2164 
2165  /* ---------------------- Monomial ----------------------- */
2166 
2167 
2168 
2169  template <int dim>
2171  const unsigned int n_components)
2172  : Function<dim>(n_components)
2173  , exponents(exponents)
2174  {}
2175 
2176 
2177 
2178  template <int dim>
2179  double
2180  Monomial<dim>::value(const Point<dim> &p, const unsigned int component) const
2181  {
2182  (void)component;
2183  AssertIndexRange(component, this->n_components);
2184 
2185  double prod = 1;
2186  for (unsigned int s = 0; s < dim; ++s)
2187  {
2188  if (p[s] < 0)
2189  Assert(std::floor(exponents[s]) == exponents[s],
2190  ExcMessage("Exponentiation of a negative base number with "
2191  "a real exponent can't be performed."));
2192  prod *= std::pow(p[s], exponents[s]);
2193  }
2194  return prod;
2195  }
2196 
2197 
2198 
2199  template <int dim>
2200  void
2202  {
2203  Assert(values.size() == this->n_components,
2204  ExcDimensionMismatch(values.size(), this->n_components));
2205 
2206  for (unsigned int i = 0; i < values.size(); ++i)
2207  values(i) = Monomial<dim>::value(p, i);
2208  }
2209 
2210 
2211 
2212  template <int dim>
2215  const unsigned int component) const
2216  {
2217  (void)component;
2218  AssertIndexRange(component, 1);
2219 
2220  Tensor<1, dim> r;
2221  for (unsigned int d = 0; d < dim; ++d)
2222  {
2223  double prod = 1;
2224  for (unsigned int s = 0; s < dim; ++s)
2225  {
2226  if ((s == d) && (exponents[s] == 0) && (p[s] == 0))
2227  {
2228  prod = 0;
2229  break;
2230  }
2231  else
2232  {
2233  if (p[s] < 0)
2234  Assert(std::floor(exponents[s]) == exponents[s],
2235  ExcMessage(
2236  "Exponentiation of a negative base number with "
2237  "a real exponent can't be performed."));
2238  prod *=
2239  (s == d ? exponents[s] * std::pow(p[s], exponents[s] - 1) :
2240  std::pow(p[s], exponents[s]));
2241  }
2242  }
2243  r[d] = prod;
2244  }
2245 
2246  return r;
2247  }
2248 
2249 
2250 
2251  template <int dim>
2252  void
2253  Monomial<dim>::value_list(const std::vector<Point<dim>> &points,
2254  std::vector<double> & values,
2255  const unsigned int component) const
2256  {
2257  Assert(values.size() == points.size(),
2258  ExcDimensionMismatch(values.size(), points.size()));
2259 
2260  for (unsigned int i = 0; i < points.size(); ++i)
2261  values[i] = Monomial<dim>::value(points[i], component);
2262  }
2263 
2264 
2265  template <int dim>
2266  Bessel1<dim>::Bessel1(const unsigned int order,
2267  const double wave_number,
2268  const Point<dim> center)
2269  : order(order)
2270  , wave_number(wave_number)
2271  , center(center)
2272  {
2273  Assert(wave_number >= 0., ExcMessage("wave_number must be nonnegative!"));
2274  }
2275 
2276  template <int dim>
2277  double
2278  Bessel1<dim>::value(const Point<dim> &p, const unsigned int) const
2279  {
2280  Assert(dim == 2, ExcNotImplemented());
2281  const double r = p.distance(center);
2282  return std_cxx17::cyl_bessel_j(order, r * wave_number);
2283  }
2284 
2285 
2286  template <int dim>
2287  void
2288  Bessel1<dim>::value_list(const std::vector<Point<dim>> &points,
2289  std::vector<double> & values,
2290  const unsigned int) const
2291  {
2292  Assert(dim == 2, ExcNotImplemented());
2293  AssertDimension(points.size(), values.size());
2294  for (unsigned int k = 0; k < points.size(); ++k)
2295  {
2296  const double r = points[k].distance(center);
2297  values[k] = std_cxx17::cyl_bessel_j(order, r * wave_number);
2298  }
2299  }
2300 
2301 
2302  template <int dim>
2304  Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
2305  {
2306  Assert(dim == 2, ExcNotImplemented());
2307  const double r = p.distance(center);
2308  const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
2309  const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
2310 
2311  const double dJn =
2312  (order == 0) ?
2313  (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2314  (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2315  std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2316  Tensor<1, dim> result;
2317  result[0] = wave_number * co * dJn;
2318  result[1] = wave_number * si * dJn;
2319  return result;
2320  }
2321 
2322 
2323 
2324  template <int dim>
2325  void
2326  Bessel1<dim>::gradient_list(const std::vector<Point<dim>> &points,
2327  std::vector<Tensor<1, dim>> & gradients,
2328  const unsigned int) const
2329  {
2330  Assert(dim == 2, ExcNotImplemented());
2331  AssertDimension(points.size(), gradients.size());
2332  for (unsigned int k = 0; k < points.size(); ++k)
2333  {
2334  const Point<dim> &p = points[k];
2335  const double r = p.distance(center);
2336  const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
2337  const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
2338 
2339  const double dJn =
2340  (order == 0) ?
2341  (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2342  (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2343  std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2344  Tensor<1, dim> &result = gradients[k];
2345  result[0] = wave_number * co * dJn;
2346  result[1] = wave_number * si * dJn;
2347  }
2348  }
2349 
2350 
2351 
2352  namespace
2353  {
2354  // interpolate a data value from a table where ix denotes
2355  // the (lower) left endpoint of the interval to interpolate
2356  // in, and p_unit denotes the point in unit coordinates to do so.
2357  double
2358  interpolate(const Table<1, double> &data_values,
2359  const TableIndices<1> & ix,
2360  const Point<1> & xi)
2361  {
2362  return ((1 - xi[0]) * data_values[ix[0]] +
2363  xi[0] * data_values[ix[0] + 1]);
2364  }
2365 
2366  double
2367  interpolate(const Table<2, double> &data_values,
2368  const TableIndices<2> & ix,
2369  const Point<2> & p_unit)
2370  {
2371  return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
2372  p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
2373  (1 - p_unit[1]) +
2374  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
2375  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
2376  p_unit[1]);
2377  }
2378 
2379  double
2380  interpolate(const Table<3, double> &data_values,
2381  const TableIndices<3> & ix,
2382  const Point<3> & p_unit)
2383  {
2384  return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
2385  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
2386  (1 - p_unit[1]) +
2387  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
2388  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
2389  p_unit[1]) *
2390  (1 - p_unit[2]) +
2391  (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
2392  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
2393  (1 - p_unit[1]) +
2394  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
2395  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
2396  p_unit[1]) *
2397  p_unit[2]);
2398  }
2399 
2400 
2401  // Interpolate the gradient of a data value from a table where ix
2402  // denotes the lower left endpoint of the interval to interpolate
2403  // in, p_unit denotes the point in unit coordinates, and dx
2404  // denotes the width of the interval in each dimension.
2405  Tensor<1, 1>
2406  gradient_interpolate(const Table<1, double> &data_values,
2407  const TableIndices<1> & ix,
2408  const Point<1> & p_unit,
2409  const Point<1> & dx)
2410  {
2411  (void)p_unit;
2412  Tensor<1, 1> grad;
2413  grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
2414  return grad;
2415  }
2416 
2417 
2418  Tensor<1, 2>
2419  gradient_interpolate(const Table<2, double> &data_values,
2420  const TableIndices<2> & ix,
2421  const Point<2> & p_unit,
2422  const Point<2> & dx)
2423  {
2424  Tensor<1, 2> grad;
2425  double u00 = data_values[ix[0]][ix[1]],
2426  u01 = data_values[ix[0] + 1][ix[1]],
2427  u10 = data_values[ix[0]][ix[1] + 1],
2428  u11 = data_values[ix[0] + 1][ix[1] + 1];
2429 
2430  grad[0] =
2431  ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
2432  grad[1] =
2433  ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
2434  return grad;
2435  }
2436 
2437 
2438  Tensor<1, 3>
2439  gradient_interpolate(const Table<3, double> &data_values,
2440  const TableIndices<3> & ix,
2441  const Point<3> & p_unit,
2442  const Point<3> & dx)
2443  {
2444  Tensor<1, 3> grad;
2445  double u000 = data_values[ix[0]][ix[1]][ix[2]],
2446  u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
2447  u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
2448  u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
2449  u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
2450  u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
2451  u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
2452  u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
2453 
2454  grad[0] =
2455  ((1 - p_unit[2]) *
2456  ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
2457  p_unit[2] *
2458  ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
2459  dx[0];
2460  grad[1] =
2461  ((1 - p_unit[2]) *
2462  ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
2463  p_unit[2] *
2464  ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
2465  dx[1];
2466  grad[2] =
2467  ((1 - p_unit[1]) *
2468  ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
2469  p_unit[1] *
2470  ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
2471  dx[2];
2472 
2473  return grad;
2474  }
2475  } // namespace
2476 
2477  template <int dim>
2479  const std::array<std::vector<double>, dim> &coordinate_values,
2480  const Table<dim, double> & data_values)
2481  : coordinate_values(coordinate_values)
2482  , data_values(data_values)
2483  {
2484  for (unsigned int d = 0; d < dim; ++d)
2485  {
2486  Assert(
2487  coordinate_values[d].size() >= 2,
2488  ExcMessage(
2489  "Coordinate arrays must have at least two coordinate values!"));
2490  for (unsigned int i = 0; i < coordinate_values[d].size() - 1; ++i)
2491  Assert(
2492  coordinate_values[d][i] < coordinate_values[d][i + 1],
2493  ExcMessage(
2494  "Coordinate arrays must be sorted in strictly ascending order."));
2495 
2496  Assert(data_values.size()[d] == coordinate_values[d].size(),
2497  ExcMessage(
2498  "Data and coordinate tables do not have the same size."));
2499  }
2500  }
2501 
2502 
2503 
2504  template <int dim>
2507  const Point<dim> &p) const
2508  {
2509  // find out where this data point lies, relative to the given
2510  // points. if we run all the way to the end of the range,
2511  // set the indices so that we will simply query the last of the
2512  // intervals, starting at x.size()-2 and going to x.size()-1.
2513  TableIndices<dim> ix;
2514  for (unsigned int d = 0; d < dim; ++d)
2515  {
2516  // get the index of the first element of the coordinate arrays that is
2517  // larger than p[d]
2518  ix[d] = (std::lower_bound(coordinate_values[d].begin(),
2519  coordinate_values[d].end(),
2520  p[d]) -
2521  coordinate_values[d].begin());
2522 
2523  // the one we want is the index of the coordinate to the left, however,
2524  // so decrease it by one (unless we have a point to the left of all, in
2525  // which case we stay where we are; the formulas below are made in a way
2526  // that allow us to extend the function by a constant value)
2527  //
2528  // to make this work, if we got coordinate_values[d].end(), we actually
2529  // have to consider the last box which has index size()-2
2530  if (ix[d] == coordinate_values[d].size())
2531  ix[d] = coordinate_values[d].size() - 2;
2532  else if (ix[d] > 0)
2533  --ix[d];
2534  }
2535 
2536  return ix;
2537  }
2538 
2539 
2540 
2541  template <int dim>
2542  double
2544  const Point<dim> & p,
2545  const unsigned int component) const
2546  {
2547  (void)component;
2548  Assert(
2549  component == 0,
2550  ExcMessage(
2551  "This is a scalar function object, the component can only be zero."));
2552 
2553  // find the index in the data table of the cell containing the input point
2554  const TableIndices<dim> ix = table_index_of_point(p);
2555 
2556  // now compute the relative point within the interval/rectangle/box
2557  // defined by the point coordinates found above. truncate below and
2558  // above to accommodate points that may lie outside the range
2559  Point<dim> p_unit;
2560  for (unsigned int d = 0; d < dim; ++d)
2561  p_unit[d] = std::max(std::min((p[d] - coordinate_values[d][ix[d]]) /
2562  (coordinate_values[d][ix[d] + 1] -
2563  coordinate_values[d][ix[d]]),
2564  1.),
2565  0.);
2566 
2567  return interpolate(data_values, ix, p_unit);
2568  }
2569 
2570 
2571 
2572  template <int dim>
2575  const Point<dim> & p,
2576  const unsigned int component) const
2577  {
2578  (void)component;
2579  Assert(
2580  component == 0,
2581  ExcMessage(
2582  "This is a scalar function object, the component can only be zero."));
2583 
2584  // find out where this data point lies
2585  const TableIndices<dim> ix = table_index_of_point(p);
2586 
2587  Point<dim> dx;
2588  for (unsigned int d = 0; d < dim; ++d)
2589  dx[d] = coordinate_values[d][ix[d] + 1] - coordinate_values[d][ix[d]];
2590 
2591  Point<dim> p_unit;
2592  for (unsigned int d = 0; d < dim; ++d)
2593  p_unit[d] =
2594  std::max(std::min((p[d] - coordinate_values[d][ix[d]]) / dx[d], 1.),
2595  0.0);
2596 
2597  return gradient_interpolate(data_values, ix, p_unit, dx);
2598  }
2599 
2600 
2601 
2602  template <int dim>
2604  const std::array<std::pair<double, double>, dim> &interval_endpoints,
2605  const std::array<unsigned int, dim> & n_subintervals,
2606  const Table<dim, double> & data_values)
2607  : interval_endpoints(interval_endpoints)
2608  , n_subintervals(n_subintervals)
2609  , data_values(data_values)
2610  {
2611  for (unsigned int d = 0; d < dim; ++d)
2612  {
2613  Assert(n_subintervals[d] >= 1,
2614  ExcMessage("There needs to be at least one subinterval in each "
2615  "coordinate direction."));
2617  ExcMessage("The interval in each coordinate direction needs "
2618  "to have positive size"));
2619  Assert(data_values.size()[d] == n_subintervals[d] + 1,
2620  ExcMessage("The data table does not have the correct size."));
2621  }
2622  }
2623 
2624 
2625  template <int dim>
2626  double
2628  const unsigned int component) const
2629  {
2630  (void)component;
2631  Assert(
2632  component == 0,
2633  ExcMessage(
2634  "This is a scalar function object, the component can only be zero."));
2635 
2636  // find out where this data point lies, relative to the given
2637  // subdivision points
2638  TableIndices<dim> ix;
2639  for (unsigned int d = 0; d < dim; ++d)
2640  {
2641  const double delta_x =
2642  ((interval_endpoints[d].second - interval_endpoints[d].first) /
2643  n_subintervals[d]);
2644  if (p[d] <= interval_endpoints[d].first)
2645  ix[d] = 0;
2646  else if (p[d] >= interval_endpoints[d].second - delta_x)
2647  ix[d] = n_subintervals[d] - 1;
2648  else
2649  ix[d] = static_cast<unsigned int>(
2650  (p[d] - interval_endpoints[d].first) / delta_x);
2651  }
2652 
2653  // now compute the relative point within the interval/rectangle/box
2654  // defined by the point coordinates found above. truncate below and
2655  // above to accommodate points that may lie outside the range
2656  Point<dim> p_unit;
2657  for (unsigned int d = 0; d < dim; ++d)
2658  {
2659  const double delta_x =
2660  ((interval_endpoints[d].second - interval_endpoints[d].first) /
2661  n_subintervals[d]);
2662  p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
2663  ix[d] * delta_x) /
2664  delta_x,
2665  1.),
2666  0.);
2667  }
2668 
2669  return interpolate(data_values, ix, p_unit);
2670  }
2671 
2672 
2673 
2674  template <int dim>
2677  const unsigned int component) const
2678  {
2679  (void)component;
2680  Assert(
2681  component == 0,
2682  ExcMessage(
2683  "This is a scalar function object, the component can only be zero."));
2684 
2685  // find out where this data point lies, relative to the given
2686  // subdivision points
2687  TableIndices<dim> ix;
2688  for (unsigned int d = 0; d < dim; ++d)
2689  {
2690  const double delta_x = ((this->interval_endpoints[d].second -
2691  this->interval_endpoints[d].first) /
2692  this->n_subintervals[d]);
2693  if (p[d] <= this->interval_endpoints[d].first)
2694  ix[d] = 0;
2695  else if (p[d] >= this->interval_endpoints[d].second - delta_x)
2696  ix[d] = this->n_subintervals[d] - 1;
2697  else
2698  ix[d] = static_cast<unsigned int>(
2699  (p[d] - this->interval_endpoints[d].first) / delta_x);
2700  }
2701 
2702  // now compute the relative point within the interval/rectangle/box
2703  // defined by the point coordinates found above. truncate below and
2704  // above to accommodate points that may lie outside the range
2705  Point<dim> p_unit;
2706  Point<dim> delta_x;
2707  for (unsigned int d = 0; d < dim; ++d)
2708  {
2709  delta_x[d] = ((this->interval_endpoints[d].second -
2710  this->interval_endpoints[d].first) /
2711  this->n_subintervals[d]);
2712  p_unit[d] =
2713  std::max(std::min((p[d] - this->interval_endpoints[d].first -
2714  ix[d] * delta_x[d]) /
2715  delta_x[d],
2716  1.),
2717  0.);
2718  }
2719 
2720  return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
2721  }
2722 
2723 
2724 
2725  /* ---------------------- Polynomial ----------------------- */
2726 
2727 
2728 
2729  template <int dim>
2731  const std::vector<double> &coefficients)
2732  : Function<dim>(1)
2733  , exponents(exponents)
2734  , coefficients(coefficients)
2735  {
2736  Assert(exponents.n_rows() == coefficients.size(),
2737  ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
2738  Assert(exponents.n_cols() == dim,
2739  ExcDimensionMismatch(exponents.n_cols(), dim));
2740  }
2741 
2742 
2743 
2744  template <int dim>
2745  double
2747  const unsigned int component) const
2748  {
2749  (void)component;
2750  AssertIndexRange(component, 1);
2751 
2752  double sum = 0;
2753  for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2754  {
2755  double prod = 1;
2756  for (unsigned int s = 0; s < dim; ++s)
2757  {
2758  if (p[s] < 0)
2759  Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
2760  ExcMessage("Exponentiation of a negative base number with "
2761  "a real exponent can't be performed."));
2762  prod *= std::pow(p[s], exponents[monom][s]);
2763  }
2764  sum += coefficients[monom] * prod;
2765  }
2766  return sum;
2767  }
2768 
2769 
2770 
2771  template <int dim>
2772  void
2773  Polynomial<dim>::value_list(const std::vector<Point<dim>> &points,
2774  std::vector<double> & values,
2775  const unsigned int component) const
2776  {
2777  Assert(values.size() == points.size(),
2778  ExcDimensionMismatch(values.size(), points.size()));
2779 
2780  for (unsigned int i = 0; i < points.size(); ++i)
2781  values[i] = Polynomial<dim>::value(points[i], component);
2782  }
2783 
2784 
2785 
2786  template <int dim>
2789  const unsigned int component) const
2790  {
2791  (void)component;
2792  AssertIndexRange(component, 1);
2793 
2794  Tensor<1, dim> r;
2795 
2796  for (unsigned int d = 0; d < dim; ++d)
2797  {
2798  double sum = 0;
2799 
2800  for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2801  {
2802  double prod = 1;
2803  for (unsigned int s = 0; s < dim; ++s)
2804  {
2805  if ((s == d) && (exponents[monom][s] == 0) && (p[s] == 0))
2806  {
2807  prod = 0;
2808  break;
2809  }
2810  else
2811  {
2812  if (p[s] < 0)
2813  Assert(std::floor(exponents[monom][s]) ==
2814  exponents[monom][s],
2815  ExcMessage(
2816  "Exponentiation of a negative base number with "
2817  "a real exponent can't be performed."));
2818  prod *=
2819  (s == d ? exponents[monom][s] *
2820  std::pow(p[s], exponents[monom][s] - 1) :
2821  std::pow(p[s], exponents[monom][s]));
2822  }
2823  }
2824  sum += coefficients[monom] * prod;
2825  }
2826  r[d] = sum;
2827  }
2828  return r;
2829  }
2830 
2831 
2832  // explicit instantiations
2833  template class SquareFunction<1>;
2834  template class SquareFunction<2>;
2835  template class SquareFunction<3>;
2836  template class Q1WedgeFunction<1>;
2837  template class Q1WedgeFunction<2>;
2838  template class Q1WedgeFunction<3>;
2839  template class PillowFunction<1>;
2840  template class PillowFunction<2>;
2841  template class PillowFunction<3>;
2842  template class CosineFunction<1>;
2843  template class CosineFunction<2>;
2844  template class CosineFunction<3>;
2845  template class CosineGradFunction<1>;
2846  template class CosineGradFunction<2>;
2847  template class CosineGradFunction<3>;
2848  template class ExpFunction<1>;
2849  template class ExpFunction<2>;
2850  template class ExpFunction<3>;
2851  template class JumpFunction<1>;
2852  template class JumpFunction<2>;
2853  template class JumpFunction<3>;
2854  template class FourierCosineFunction<1>;
2855  template class FourierCosineFunction<2>;
2856  template class FourierCosineFunction<3>;
2857  template class FourierSineFunction<1>;
2858  template class FourierSineFunction<2>;
2859  template class FourierSineFunction<3>;
2860  template class FourierCosineSum<1>;
2861  template class FourierCosineSum<2>;
2862  template class FourierCosineSum<3>;
2863  template class FourierSineSum<1>;
2864  template class FourierSineSum<2>;
2865  template class FourierSineSum<3>;
2866  template class SlitSingularityFunction<2>;
2867  template class SlitSingularityFunction<3>;
2868  template class Monomial<1>;
2869  template class Monomial<2>;
2870  template class Monomial<3>;
2871  template class Bessel1<1>;
2872  template class Bessel1<2>;
2873  template class Bessel1<3>;
2874  template class InterpolatedTensorProductGridData<1>;
2875  template class InterpolatedTensorProductGridData<2>;
2876  template class InterpolatedTensorProductGridData<3>;
2877  template class InterpolatedUniformGridData<1>;
2878  template class InterpolatedUniformGridData<2>;
2879  template class InterpolatedUniformGridData<3>;
2880  template class Polynomial<1>;
2881  template class Polynomial<2>;
2882  template class Polynomial<3>;
2883 } // namespace Functions
2884 
Functions::Monomial
Definition: function_lib.h:1341
Functions::CosineGradFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:748
Functions::CosineFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:556
Functions::InterpolatedTensorProductGridData::coordinate_values
const std::array< std::vector< double >, dim > coordinate_values
Definition: function_lib.h:1484
Functions::SquareFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:69
Functions::CosineGradFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:868
Functions
Definition: flow_function.h:28
Differentiation::SD::floor
Expression floor(const Expression &x)
Definition: symengine_math.cc:294
Functions::CosineGradFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
Definition: function_lib.cc:955
Functions::LSingularityFunction::gradient
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1263
Functions::LSingularityGradFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< 2 >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:1385
Functions::FourierCosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1956
Functions::Monomial::Monomial
Monomial(const Tensor< 1, dim > &exponents, const unsigned int n_components=1)
Definition: function_lib.cc:2170
Functions::Monomial::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:2253
Functions::FourierSineFunction::FourierSineFunction
FourierSineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition: function_lib.cc:1972
Functions::InterpolatedUniformGridData
Definition: function_lib.h:1530
Functions::CosineFunction
Definition: function_lib.h:223
Functions::CosineFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:468
Functions::FourierCosineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2112
Functions::SlitSingularityFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1527
Functions::CosineFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:510
Functions::JumpFunction::memory_consumption
std::size_t memory_consumption() const
Definition: function_lib.cc:1909
TableIndices
Definition: table_indices.h:45
Functions::SquareFunction::vector_gradient
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim >> &gradient) const override
Definition: function_lib.cc:100
Functions::InterpolatedUniformGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2676
Functions::InterpolatedTensorProductGridData::table_index_of_point
TableIndices< dim > table_index_of_point(const Point< dim > &p) const
Definition: function_lib.cc:2506
SymmetricTensor< 2, dim >
Functions::FourierSineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1994
Differentiation::SD::atan
Expression atan(const Expression &x)
Definition: symengine_math.cc:147
Functions::FourierCosineSum::FourierCosineSum
FourierCosineSum(const std::vector< Point< dim >> &fourier_coefficients, const std::vector< double > &weights)
Definition: function_lib.cc:2096
Functions::SquareFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:52
Functions::SquareFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:34
Functions::ExpFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1037
Functions::Monomial::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:2201
Functions::InterpolatedTensorProductGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2543
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Functions::Q1WedgeFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:129
Functions::LSingularityFunction::value_list
virtual void value_list(const std::vector< Point< 2 >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1189
Functions::Q1WedgeFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< dim >> &, std::vector< std::vector< Tensor< 1, dim >>> &) const override
Definition: function_lib.cc:232
Functions::InterpolatedTensorProductGridData
Definition: function_lib.h:1426
Functions::InterpolatedTensorProductGridData::data_values
const Table< dim, double > data_values
Definition: function_lib.h:1489
DataOutBase::dx
@ dx
Definition: data_out_base.h:1562
Functions::FourierSineSum::weights
const std::vector< double > weights
Definition: function_lib.h:852
Functions::Bessel1::wave_number
double wave_number
Definition: function_bessel.h:67
Functions::JumpFunction::sine
double sine
Definition: function_lib.h:684
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
Functions::SlitHyperSingularityFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< 2 >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:1669
Functions::FourierSineSum::FourierSineSum
FourierSineSum(const std::vector< Point< dim >> &fourier_coefficients, const std::vector< double > &weights)
Definition: function_lib.cc:2022
Functions::Bessel1::value
virtual double value(const Point< dim > &points, const unsigned int component=0) const override
Definition: function_lib.cc:2278
Functions::JumpFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1827
Functions::FourierCosineFunction::FourierCosineFunction
FourierCosineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition: function_lib.cc:1922
Functions::SlitHyperSingularityFunction::value_list
virtual void value_list(const std::vector< Point< 2 >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1648
Functions::CosineFunction::CosineFunction
CosineFunction(const unsigned int n_components=1)
Definition: function_lib.cc:425
Functions::CosineFunction::hessian
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:598
Vector::size
size_type size() const
Functions::FourierCosineFunction
Definition: function_lib.h:707
Functions::FourierSineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2074
Functions::CosineGradFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:859
Functions::Q1WedgeFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:157
Functions::ExpFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1084
Functions::Polynomial::coefficients
const std::vector< double > coefficients
Definition: function_lib.h:1658
second
Point< 2 > second
Definition: grid_out.cc:4353
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Functions::PillowFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:333
Functions::CosineGradFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component) const override
Definition: function_lib.cc:908
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Functions::FourierSineSum
Definition: function_lib.h:813
Functions::FourierCosineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2130
numbers::PI_2
static constexpr double PI_2
Definition: numbers.h:242
Functions::PillowFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:280
cmath.h
Functions::LSingularityGradFunction::value_list
virtual void value_list(const std::vector< Point< 2 >> &points, std::vector< double > &values, const unsigned int component) const override
Definition: function_lib.cc:1360
Functions::Bessel1::Bessel1
Bessel1(const unsigned int order, const double wave_number, const Point< dim > center=Point< dim >())
Definition: function_lib.cc:2266
Table
Definition: table.h:699
Functions::SlitSingularityFunction
Definition: function_lib.h:500
Functions::SlitSingularityFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1551
Functions::ExpFunction
Definition: function_lib.h:339
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1607
Functions::SlitHyperSingularityFunction::value
virtual double value(const Point< 2 > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1634
Functions::FourierCosineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition: function_lib.h:906
Functions::FourierCosineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2148
Differentiation::SD::atan2
Expression atan2(const Expression &y, const Expression &x)
Definition: symengine_math.cc:154
Functions::JumpFunction::JumpFunction
JumpFunction(const Point< dim > &direction, const double steepness)
Definition: function_lib.cc:1793
function_bessel.h
Functions::SlitSingularityFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1464
Functions::LSingularityGradFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< 2 >> &, std::vector< std::vector< Tensor< 1, 2 >>> &) const override
Definition: function_lib.cc:1453
function_lib.h
tensor.h
Functions::Monomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2214
Functions::SquareFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:92
Functions::LSingularityFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< 2 >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:1215
Functions::JumpFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:1889
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
Functions::JumpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1844
Functions::SlitSingularityFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:1502
Functions::CosineFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:454
Tensor< 1, dim >
Functions::PillowFunction
Definition: function_lib.h:155
Functions::LSingularityGradFunction::laplacian
virtual double laplacian(const Point< 2 > &p, const unsigned int component) const override
Definition: function_lib.cc:1412
Functions::Polynomial
Definition: function_lib.h:1612
Utilities::lower_bound
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1102
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
Functions::SlitHyperSingularityFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< 2 >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1701
Functions::FourierCosineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1932
Functions::FourierSineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2006
Functions::ExpFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1113
Functions::Bessel1::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2304
Functions::SquareFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:42
Functions::Polynomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2788
Functions::Polynomial::Polynomial
Polynomial(const Table< 2, double > &exponents, const std::vector< double > &coefficients)
Definition: function_lib.cc:2730
Functions::FourierSineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2038
Functions::LSingularityGradFunction::LSingularityGradFunction
LSingularityGradFunction()
Definition: function_lib.cc:1336
Functions::CosineFunction::hessian_list
virtual void hessian_list(const std::vector< Point< dim >> &points, std::vector< SymmetricTensor< 2, dim >> &hessians, const unsigned int component=0) const override
Definition: function_lib.cc:652
Functions::FourierSineFunction
Definition: function_lib.h:762
Functions::CosineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:433
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
Functions::ExpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1066
Functions::JumpFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1817
Functions::LSingularityGradFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< 2 >> &points, std::vector< double > &values, const unsigned int component) const override
Definition: function_lib.cc:1420
Functions::FourierSineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1982
Functions::InterpolatedTensorProductGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2574
TableBase::size
size_type size(const unsigned int i) const
Functions::JumpFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1875
Functions::LSingularityGradFunction::gradient
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component) const override
Definition: function_lib.cc:1434
Functions::LSingularityFunction::laplacian
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1243
Functions::Q1WedgeFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:214
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Functions::LSingularityGradFunction::gradient_list
virtual void gradient_list(const std::vector< Point< 2 >> &points, std::vector< Tensor< 1, 2 >> &gradients, const unsigned int component) const override
Definition: function_lib.cc:1443
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Functions::InterpolatedUniformGridData::n_subintervals
const std::array< unsigned int, dim > n_subintervals
Definition: function_lib.h:1589
Functions::CosineGradFunction
Definition: function_lib.h:292
Functions::Bessel1::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:2326
Functions::LSingularityGradFunction::value
virtual double value(const Point< 2 > &p, const unsigned int component) const override
Definition: function_lib.cc:1342
Functions::Bessel1::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:2288
Functions::InterpolatedUniformGridData::data_values
const Table< dim, double > data_values
Definition: function_lib.h:1594
Functions::InterpolatedUniformGridData::interval_endpoints
const std::array< std::pair< double, double >, dim > interval_endpoints
Definition: function_lib.h:1584
Functions::PillowFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:364
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Functions::SlitHyperSingularityFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< 2 >> &, std::vector< std::vector< Tensor< 1, 2 >>> &) const override
Definition: function_lib.cc:1763
Functions::PillowFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:313
vector.h
Functions::Q1WedgeFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:176
Functions::Q1WedgeFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:139
Functions::Polynomial::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2746
Functions::Q1WedgeFunction
Definition: function_lib.h:92
Functions::Q1WedgeFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:201
Functions::SlitSingularityFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1536
Functions::Q1WedgeFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:185
Functions::Monomial::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2180
Functions::SlitHyperSingularityFunction::gradient
virtual Tensor< 1, 2 > gradient(const Point< 2 > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1715
Functions::Polynomial::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:2773
Functions::JumpFunction::cosine
double cosine
Definition: function_lib.h:689
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
StandardExceptions::ExcZero
static ::ExceptionBase & ExcZero()
Functions::CosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:523
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
Functions::CosineGradFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component) const override
Definition: function_lib.cc:782
Functions::PillowFunction::PillowFunction
PillowFunction(const double offset=0.)
Definition: function_lib.cc:254
Functions::ExpFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:1138
Functions::CosineGradFunction::value
virtual double value(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:722
Functions::FourierSineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2056
Functions::JumpFunction::angle
double angle
Definition: function_lib.h:679
Point::square
numbers::NumberTraits< Number >::real_type square() const
Functions::SlitSingularityFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1479
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
Functions::JumpFunction::direction
const Point< dim > direction
Definition: function_lib.h:669
Function
Definition: function.h:151
Functions::PillowFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:261
Functions::SlitSingularityFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< dim >> &, std::vector< std::vector< Tensor< 1, dim >>> &) const override
Definition: function_lib.cc:1601
Functions::JumpFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1854
Functions::SlitSingularityFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:1572
Functions::InterpolatedUniformGridData::InterpolatedUniformGridData
InterpolatedUniformGridData(const std::array< std::pair< double, double >, dim > &interval_endpoints, const std::array< unsigned int, dim > &n_subintervals, const Table< dim, double > &data_values)
Definition: function_lib.cc:2603
Functions::SquareFunction
Definition: function_lib.h:52
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Functions::InterpolatedTensorProductGridData::InterpolatedTensorProductGridData
InterpolatedTensorProductGridData(const std::array< std::vector< double >, dim > &coordinate_values, const Table< dim, double > &data_values)
Definition: function_lib.cc:2478
Functions::InterpolatedUniformGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2627
Functions::SquareFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:77
first
Point< 2 > first
Definition: grid_out.cc:4352
Functions::FourierSineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition: function_lib.h:851
Functions::CosineGradFunction::CosineGradFunction
CosineGradFunction()
Definition: function_lib.cc:715
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Functions::PillowFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:389
Vector< double >
Functions::FourierCosineSum
Definition: function_lib.h:868
Functions::FourierCosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1944
Functions::SlitHyperSingularityFunction::laplacian
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1693
std_cxx17::cyl_bessel_j
double(&) cyl_bessel_j(double, double)
Definition: cmath.h:30
Functions::LSingularityFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< 2 >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:1250
Functions::FourierCosineSum::weights
const std::vector< double > weights
Definition: function_lib.h:907
Functions::SlitHyperSingularityFunction::gradient_list
virtual void gradient_list(const std::vector< Point< 2 >> &points, std::vector< Tensor< 1, 2 >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:1736
center
Point< 3 > center
Definition: data_out_base.cc:169
numbers.h
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
Functions::LSingularityFunction::vector_gradient_list
virtual void vector_gradient_list(const std::vector< Point< 2 >> &, std::vector< std::vector< Tensor< 1, 2 >>> &) const override
Definition: function_lib.cc:1308
Functions::CosineGradFunction::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
Definition: function_lib.cc:818
Functions::LSingularityFunction::gradient_list
virtual void gradient_list(const std::vector< Point< 2 >> &points, std::vector< Tensor< 1, 2 >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:1282
Functions::CosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:486
Functions::LSingularityFunction::value
virtual double value(const Point< 2 > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1173
Functions::Bessel1
Definition: function_bessel.h:37
Functions::Polynomial::exponents
const Table< 2, double > exponents
Definition: function_lib.h:1653
Functions::SquareFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:112
Functions::JumpFunction
Definition: function_lib.h:603
point.h
Functions::ExpFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1019