Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
function_lib.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/function_bessel.h>
17 #include <deal.II/base/function_lib.h>
18 #include <deal.II/base/numbers.h>
19 #include <deal.II/base/point.h>
20 #include <deal.II/base/std_cxx17/cmath.h>
21 #include <deal.II/base/tensor.h>
22 
23 #include <deal.II/lac/vector.h>
24 
25 #include <cmath>
26 
27 DEAL_II_NAMESPACE_OPEN
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
42  SquareFunction<dim>::vector_value(const Point<dim> &p,
43  Vector<double> & values) const
44  {
45  AssertDimension(values.size(), 1);
46  values(0) = p.square();
47  }
48 
49 
50  template <int dim>
51  void
52  SquareFunction<dim>::value_list(const std::vector<Point<dim>> &points,
53  std::vector<double> & values,
54  const unsigned int) const
55  {
56  Assert(values.size() == points.size(),
57  ExcDimensionMismatch(values.size(), points.size()));
58 
59  for (unsigned int i = 0; i < points.size(); ++i)
60  {
61  const Point<dim> &p = points[i];
62  values[i] = p.square();
63  }
64  }
65 
66 
67  template <int dim>
68  double
69  SquareFunction<dim>::laplacian(const Point<dim> &, const unsigned int) const
70  {
71  return 2 * dim;
72  }
73 
74 
75  template <int dim>
76  void
77  SquareFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
78  std::vector<double> & values,
79  const unsigned int) const
80  {
81  Assert(values.size() == points.size(),
82  ExcDimensionMismatch(values.size(), points.size()));
83 
84  for (unsigned int i = 0; i < points.size(); ++i)
85  values[i] = 2 * dim;
86  }
87 
88 
89 
90  template <int dim>
91  Tensor<1, dim>
92  SquareFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
93  {
94  return p * 2.;
95  }
96 
97 
98  template <int dim>
99  void
100  SquareFunction<dim>::vector_gradient(
101  const Point<dim> & p,
102  std::vector<Tensor<1, dim>> &values) const
103  {
104  AssertDimension(values.size(), 1);
105  values[0] = p * 2.;
106  }
107 
108 
109 
110  template <int dim>
111  void
112  SquareFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
113  std::vector<Tensor<1, dim>> & gradients,
114  const unsigned int) const
115  {
116  Assert(gradients.size() == points.size(),
117  ExcDimensionMismatch(gradients.size(), points.size()));
118 
119  for (unsigned int i = 0; i < points.size(); ++i)
120  gradients[i] = static_cast<Tensor<1, dim>>(points[i]) * 2;
121  }
122 
123 
125 
126 
127  template <int dim>
128  double
129  Q1WedgeFunction<dim>::value(const Point<dim> &p, const unsigned int) const
130  {
131  Assert(dim >= 2, ExcInternalError());
132  return p(0) * p(1);
133  }
134 
135 
136 
137  template <int dim>
138  void
139  Q1WedgeFunction<dim>::value_list(const std::vector<Point<dim>> &points,
140  std::vector<double> & values,
141  const unsigned int) const
142  {
143  Assert(dim >= 2, ExcInternalError());
144  Assert(values.size() == points.size(),
145  ExcDimensionMismatch(values.size(), points.size()));
146 
147  for (unsigned int i = 0; i < points.size(); ++i)
148  {
149  const Point<dim> &p = points[i];
150  values[i] = p(0) * p(1);
151  }
152  }
153 
154 
155  template <int dim>
156  void
157  Q1WedgeFunction<dim>::vector_value_list(
158  const std::vector<Point<dim>> &points,
159  std::vector<Vector<double>> & values) const
160  {
161  Assert(dim >= 2, ExcInternalError());
162  Assert(values.size() == points.size(),
163  ExcDimensionMismatch(values.size(), points.size()));
164  Assert(values[0].size() == 1, ExcDimensionMismatch(values[0].size(), 1));
165 
166  for (unsigned int i = 0; i < points.size(); ++i)
167  {
168  const Point<dim> &p = points[i];
169  values[i](0) = p(0) * p(1);
170  }
171  }
172 
173 
174  template <int dim>
175  double
176  Q1WedgeFunction<dim>::laplacian(const Point<dim> &, const unsigned int) const
177  {
178  Assert(dim >= 2, ExcInternalError());
179  return 0.;
180  }
181 
182 
183  template <int dim>
184  void
185  Q1WedgeFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
186  std::vector<double> & values,
187  const unsigned int) const
188  {
189  Assert(dim >= 2, ExcInternalError());
190  Assert(values.size() == points.size(),
191  ExcDimensionMismatch(values.size(), points.size()));
192 
193  for (unsigned int i = 0; i < points.size(); ++i)
194  values[i] = 0.;
195  }
196 
197 
198 
199  template <int dim>
200  Tensor<1, dim>
201  Q1WedgeFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
202  {
203  Assert(dim >= 2, ExcInternalError());
204  Tensor<1, dim> erg;
205  erg[0] = p(1);
206  erg[1] = p(0);
207  return erg;
208  }
209 
210 
211 
212  template <int dim>
213  void
214  Q1WedgeFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
215  std::vector<Tensor<1, dim>> & gradients,
216  const unsigned int) const
217  {
218  Assert(dim >= 2, ExcInternalError());
219  Assert(gradients.size() == points.size(),
220  ExcDimensionMismatch(gradients.size(), points.size()));
221 
222  for (unsigned int i = 0; i < points.size(); ++i)
223  {
224  gradients[i][0] = points[i](1);
225  gradients[i][1] = points[i](0);
226  }
227  }
228 
229 
230  template <int dim>
231  void
232  Q1WedgeFunction<dim>::vector_gradient_list(
233  const std::vector<Point<dim>> & points,
234  std::vector<std::vector<Tensor<1, dim>>> &gradients) const
235  {
236  Assert(dim >= 2, ExcInternalError());
237  Assert(gradients.size() == points.size(),
238  ExcDimensionMismatch(gradients.size(), points.size()));
239  Assert(gradients[0].size() == 1,
240  ExcDimensionMismatch(gradients[0].size(), 1));
241 
242  for (unsigned int i = 0; i < points.size(); ++i)
243  {
244  gradients[i][0][0] = points[i](1);
245  gradients[i][0][1] = points[i](0);
246  }
247  }
248 
249 
251 
252 
253  template <int dim>
254  PillowFunction<dim>::PillowFunction(const double offset)
255  : offset(offset)
256  {}
257 
258 
259  template <int dim>
260  double
261  PillowFunction<dim>::value(const Point<dim> &p, const unsigned int) const
262  {
263  switch (dim)
264  {
265  case 1:
266  return 1. - p(0) * p(0) + offset;
267  case 2:
268  return (1. - p(0) * p(0)) * (1. - p(1) * p(1)) + offset;
269  case 3:
270  return (1. - p(0) * p(0)) * (1. - p(1) * p(1)) * (1. - p(2) * p(2)) +
271  offset;
272  default:
273  Assert(false, ExcNotImplemented());
274  }
275  return 0.;
276  }
277 
278  template <int dim>
279  void
280  PillowFunction<dim>::value_list(const std::vector<Point<dim>> &points,
281  std::vector<double> & values,
282  const unsigned int) const
283  {
284  Assert(values.size() == points.size(),
285  ExcDimensionMismatch(values.size(), points.size()));
286 
287  for (unsigned int i = 0; i < points.size(); ++i)
288  {
289  const Point<dim> &p = points[i];
290  switch (dim)
291  {
292  case 1:
293  values[i] = 1. - p(0) * p(0) + offset;
294  break;
295  case 2:
296  values[i] = (1. - p(0) * p(0)) * (1. - p(1) * p(1)) + offset;
297  break;
298  case 3:
299  values[i] =
300  (1. - p(0) * p(0)) * (1. - p(1) * p(1)) * (1. - p(2) * p(2)) +
301  offset;
302  break;
303  default:
304  Assert(false, ExcNotImplemented());
305  }
306  }
307  }
308 
309 
310 
311  template <int dim>
312  double
313  PillowFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
314  {
315  switch (dim)
316  {
317  case 1:
318  return -2.;
319  case 2:
320  return -2. * ((1. - p(0) * p(0)) + (1. - p(1) * p(1)));
321  case 3:
322  return -2. * ((1. - p(0) * p(0)) * (1. - p(1) * p(1)) +
323  (1. - p(1) * p(1)) * (1. - p(2) * p(2)) +
324  (1. - p(2) * p(2)) * (1. - p(0) * p(0)));
325  default:
326  Assert(false, ExcNotImplemented());
327  }
328  return 0.;
329  }
330 
331  template <int dim>
332  void
333  PillowFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
334  std::vector<double> & values,
335  const unsigned int) const
336  {
337  Assert(values.size() == points.size(),
338  ExcDimensionMismatch(values.size(), points.size()));
339 
340  for (unsigned int i = 0; i < points.size(); ++i)
341  {
342  const Point<dim> &p = points[i];
343  switch (dim)
344  {
345  case 1:
346  values[i] = -2.;
347  break;
348  case 2:
349  values[i] = -2. * ((1. - p(0) * p(0)) + (1. - p(1) * p(1)));
350  break;
351  case 3:
352  values[i] = -2. * ((1. - p(0) * p(0)) * (1. - p(1) * p(1)) +
353  (1. - p(1) * p(1)) * (1. - p(2) * p(2)) +
354  (1. - p(2) * p(2)) * (1. - p(0) * p(0)));
355  break;
356  default:
357  Assert(false, ExcNotImplemented());
358  }
359  }
360  }
361 
362  template <int dim>
363  Tensor<1, dim>
364  PillowFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
365  {
366  Tensor<1, dim> result;
367  switch (dim)
368  {
369  case 1:
370  result[0] = -2. * p(0);
371  break;
372  case 2:
373  result[0] = -2. * p(0) * (1. - p(1) * p(1));
374  result[1] = -2. * p(1) * (1. - p(0) * p(0));
375  break;
376  case 3:
377  result[0] = -2. * p(0) * (1. - p(1) * p(1)) * (1. - p(2) * p(2));
378  result[1] = -2. * p(1) * (1. - p(0) * p(0)) * (1. - p(2) * p(2));
379  result[2] = -2. * p(2) * (1. - p(0) * p(0)) * (1. - p(1) * p(1));
380  break;
381  default:
382  Assert(false, ExcNotImplemented());
383  }
384  return result;
385  }
386 
387  template <int dim>
388  void
389  PillowFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
390  std::vector<Tensor<1, dim>> & gradients,
391  const unsigned int) const
392  {
393  Assert(gradients.size() == points.size(),
394  ExcDimensionMismatch(gradients.size(), points.size()));
395 
396  for (unsigned int i = 0; i < points.size(); ++i)
397  {
398  const Point<dim> &p = points[i];
399  switch (dim)
400  {
401  case 1:
402  gradients[i][0] = -2. * p(0);
403  break;
404  case 2:
405  gradients[i][0] = -2. * p(0) * (1. - p(1) * p(1));
406  gradients[i][1] = -2. * p(1) * (1. - p(0) * p(0));
407  break;
408  case 3:
409  gradients[i][0] =
410  -2. * p(0) * (1. - p(1) * p(1)) * (1. - p(2) * p(2));
411  gradients[i][1] =
412  -2. * p(1) * (1. - p(0) * p(0)) * (1. - p(2) * p(2));
413  gradients[i][2] =
414  -2. * p(2) * (1. - p(0) * p(0)) * (1. - p(1) * p(1));
415  break;
416  default:
417  Assert(false, ExcNotImplemented());
418  }
419  }
420  }
421 
423 
424  template <int dim>
425  CosineFunction<dim>::CosineFunction(const unsigned int n_components)
426  : Function<dim>(n_components)
427  {}
428 
429 
430 
431  template <int dim>
432  double
433  CosineFunction<dim>::value(const Point<dim> &p, const unsigned int) const
434  {
435  switch (dim)
436  {
437  case 1:
438  return std::cos(numbers::PI_2 * p(0));
439  case 2:
440  return std::cos(numbers::PI_2 * p(0)) *
441  std::cos(numbers::PI_2 * p(1));
442  case 3:
443  return std::cos(numbers::PI_2 * p(0)) *
444  std::cos(numbers::PI_2 * p(1)) *
445  std::cos(numbers::PI_2 * p(2));
446  default:
447  Assert(false, ExcNotImplemented());
448  }
449  return 0.;
450  }
451 
452  template <int dim>
453  void
454  CosineFunction<dim>::value_list(const std::vector<Point<dim>> &points,
455  std::vector<double> & values,
456  const unsigned int) const
457  {
458  Assert(values.size() == points.size(),
459  ExcDimensionMismatch(values.size(), points.size()));
460 
461  for (unsigned int i = 0; i < points.size(); ++i)
462  values[i] = value(points[i]);
463  }
464 
465 
466  template <int dim>
467  void
468  CosineFunction<dim>::vector_value_list(
469  const std::vector<Point<dim>> &points,
470  std::vector<Vector<double>> & values) const
471  {
472  Assert(values.size() == points.size(),
473  ExcDimensionMismatch(values.size(), points.size()));
474 
475  for (unsigned int i = 0; i < points.size(); ++i)
476  {
477  const double v = value(points[i]);
478  for (unsigned int k = 0; k < values[i].size(); ++k)
479  values[i](k) = v;
480  }
481  }
482 
483 
484  template <int dim>
485  double
486  CosineFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
487  {
488  switch (dim)
489  {
490  case 1:
491  return -numbers::PI_2 * numbers::PI_2 *
492  std::cos(numbers::PI_2 * p(0));
493  case 2:
494  return -2 * numbers::PI_2 * numbers::PI_2 *
495  std::cos(numbers::PI_2 * p(0)) *
496  std::cos(numbers::PI_2 * p(1));
497  case 3:
498  return -3 * numbers::PI_2 * numbers::PI_2 *
499  std::cos(numbers::PI_2 * p(0)) *
500  std::cos(numbers::PI_2 * p(1)) *
501  std::cos(numbers::PI_2 * p(2));
502  default:
503  Assert(false, ExcNotImplemented());
504  }
505  return 0.;
506  }
507 
508  template <int dim>
509  void
510  CosineFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
511  std::vector<double> & values,
512  const unsigned int) const
513  {
514  Assert(values.size() == points.size(),
515  ExcDimensionMismatch(values.size(), points.size()));
516 
517  for (unsigned int i = 0; i < points.size(); ++i)
518  values[i] = laplacian(points[i]);
519  }
520 
521  template <int dim>
522  Tensor<1, dim>
523  CosineFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
524  {
525  Tensor<1, dim> result;
526  switch (dim)
527  {
528  case 1:
529  result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0));
530  break;
531  case 2:
532  result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
533  std::cos(numbers::PI_2 * p(1));
534  result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
535  std::sin(numbers::PI_2 * p(1));
536  break;
537  case 3:
538  result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
539  std::cos(numbers::PI_2 * p(1)) *
540  std::cos(numbers::PI_2 * p(2));
541  result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
542  std::sin(numbers::PI_2 * p(1)) *
543  std::cos(numbers::PI_2 * p(2));
544  result[2] = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
545  std::cos(numbers::PI_2 * p(1)) *
546  std::sin(numbers::PI_2 * p(2));
547  break;
548  default:
549  Assert(false, ExcNotImplemented());
550  }
551  return result;
552  }
553 
554  template <int dim>
555  void
556  CosineFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
557  std::vector<Tensor<1, dim>> & gradients,
558  const unsigned int) const
559  {
560  Assert(gradients.size() == points.size(),
561  ExcDimensionMismatch(gradients.size(), points.size()));
562 
563  for (unsigned int i = 0; i < points.size(); ++i)
564  {
565  const Point<dim> &p = points[i];
566  switch (dim)
567  {
568  case 1:
569  gradients[i][0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0));
570  break;
571  case 2:
572  gradients[i][0] = -numbers::PI_2 *
573  std::sin(numbers::PI_2 * p(0)) *
574  std::cos(numbers::PI_2 * p(1));
575  gradients[i][1] = -numbers::PI_2 *
576  std::cos(numbers::PI_2 * p(0)) *
577  std::sin(numbers::PI_2 * p(1));
578  break;
579  case 3:
580  gradients[i][0] =
581  -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
582  std::cos(numbers::PI_2 * p(1)) * std::cos(numbers::PI_2 * p(2));
583  gradients[i][1] =
584  -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
585  std::sin(numbers::PI_2 * p(1)) * std::cos(numbers::PI_2 * p(2));
586  gradients[i][2] =
587  -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
588  std::cos(numbers::PI_2 * p(1)) * std::sin(numbers::PI_2 * p(2));
589  break;
590  default:
591  Assert(false, ExcNotImplemented());
592  }
593  }
594  }
595 
596  template <int dim>
597  SymmetricTensor<2, dim>
598  CosineFunction<dim>::hessian(const Point<dim> &p, const unsigned int) const
599  {
600  const double pi2 = numbers::PI_2 * numbers::PI_2;
601 
602  SymmetricTensor<2, dim> result;
603  switch (dim)
604  {
605  case 1:
606  result[0][0] = -pi2 * std::cos(numbers::PI_2 * p(0));
607  break;
608  case 2:
609  {
610  const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) *
611  std::cos(numbers::PI_2 * p(1));
612  const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) *
613  std::sin(numbers::PI_2 * p(1));
614  result[0][0] = coco;
615  result[1][1] = coco;
616  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
617  result[0][1] = sisi;
618  }
619  break;
620  case 3:
621  {
622  const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) *
623  std::cos(numbers::PI_2 * p(1)) *
624  std::cos(numbers::PI_2 * p(2));
625  const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) *
626  std::sin(numbers::PI_2 * p(1)) *
627  std::cos(numbers::PI_2 * p(2));
628  const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) *
629  std::cos(numbers::PI_2 * p(1)) *
630  std::sin(numbers::PI_2 * p(2));
631  const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) *
632  std::sin(numbers::PI_2 * p(1)) *
633  std::sin(numbers::PI_2 * p(2));
634 
635  result[0][0] = cococo;
636  result[1][1] = cococo;
637  result[2][2] = cococo;
638  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
639  result[0][1] = sisico;
640  result[0][2] = sicosi;
641  result[1][2] = cosisi;
642  }
643  break;
644  default:
645  Assert(false, ExcNotImplemented());
646  }
647  return result;
648  }
649 
650  template <int dim>
651  void
652  CosineFunction<dim>::hessian_list(
653  const std::vector<Point<dim>> & points,
654  std::vector<SymmetricTensor<2, dim>> &hessians,
655  const unsigned int) const
656  {
657  Assert(hessians.size() == points.size(),
658  ExcDimensionMismatch(hessians.size(), points.size()));
659 
660  const double pi2 = numbers::PI_2 * numbers::PI_2;
661 
662  for (unsigned int i = 0; i < points.size(); ++i)
663  {
664  const Point<dim> &p = points[i];
665  switch (dim)
666  {
667  case 1:
668  hessians[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p(0));
669  break;
670  case 2:
671  {
672  const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) *
673  std::cos(numbers::PI_2 * p(1));
674  const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) *
675  std::sin(numbers::PI_2 * p(1));
676  hessians[i][0][0] = coco;
677  hessians[i][1][1] = coco;
678  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
679  hessians[i][0][1] = sisi;
680  }
681  break;
682  case 3:
683  {
684  const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) *
685  std::cos(numbers::PI_2 * p(1)) *
686  std::cos(numbers::PI_2 * p(2));
687  const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) *
688  std::sin(numbers::PI_2 * p(1)) *
689  std::cos(numbers::PI_2 * p(2));
690  const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) *
691  std::cos(numbers::PI_2 * p(1)) *
692  std::sin(numbers::PI_2 * p(2));
693  const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) *
694  std::sin(numbers::PI_2 * p(1)) *
695  std::sin(numbers::PI_2 * p(2));
696 
697  hessians[i][0][0] = cococo;
698  hessians[i][1][1] = cococo;
699  hessians[i][2][2] = cococo;
700  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
701  hessians[i][0][1] = sisico;
702  hessians[i][0][2] = sicosi;
703  hessians[i][1][2] = cosisi;
704  }
705  break;
706  default:
707  Assert(false, ExcNotImplemented());
708  }
709  }
710  }
711 
713 
714  template <int dim>
715  CosineGradFunction<dim>::CosineGradFunction()
716  : Function<dim>(dim)
717  {}
718 
719 
720  template <int dim>
721  double
722  CosineGradFunction<dim>::value(const Point<dim> & p,
723  const unsigned int d) const
724  {
725  AssertIndexRange(d, dim);
726  const unsigned int d1 = (d + 1) % dim;
727  const unsigned int d2 = (d + 2) % dim;
728  switch (dim)
729  {
730  case 1:
731  return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(0)));
732  case 2:
733  return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
734  std::cos(numbers::PI_2 * p(d1)));
735  case 3:
736  return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
737  std::cos(numbers::PI_2 * p(d1)) *
738  std::cos(numbers::PI_2 * p(d2)));
739  default:
740  Assert(false, ExcNotImplemented());
741  }
742  return 0.;
743  }
744 
745 
746  template <int dim>
747  void
748  CosineGradFunction<dim>::vector_value(const Point<dim> &p,
749  Vector<double> & result) const
750  {
751  AssertDimension(result.size(), dim);
752  switch (dim)
753  {
754  case 1:
755  result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0));
756  break;
757  case 2:
758  result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
759  std::cos(numbers::PI_2 * p(1));
760  result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
761  std::sin(numbers::PI_2 * p(1));
762  break;
763  case 3:
764  result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
765  std::cos(numbers::PI_2 * p(1)) *
766  std::cos(numbers::PI_2 * p(2));
767  result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
768  std::sin(numbers::PI_2 * p(1)) *
769  std::cos(numbers::PI_2 * p(2));
770  result(2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
771  std::cos(numbers::PI_2 * p(1)) *
772  std::sin(numbers::PI_2 * p(2));
773  break;
774  default:
775  Assert(false, ExcNotImplemented());
776  }
777  }
778 
779 
780  template <int dim>
781  void
782  CosineGradFunction<dim>::value_list(const std::vector<Point<dim>> &points,
783  std::vector<double> & values,
784  const unsigned int d) const
785  {
786  Assert(values.size() == points.size(),
787  ExcDimensionMismatch(values.size(), points.size()));
788  AssertIndexRange(d, dim);
789  const unsigned int d1 = (d + 1) % dim;
790  const unsigned int d2 = (d + 2) % dim;
791 
792  for (unsigned int i = 0; i < points.size(); ++i)
793  {
794  const Point<dim> &p = points[i];
795  switch (dim)
796  {
797  case 1:
798  values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d));
799  break;
800  case 2:
801  values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
802  std::cos(numbers::PI_2 * p(d1));
803  break;
804  case 3:
805  values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
806  std::cos(numbers::PI_2 * p(d1)) *
807  std::cos(numbers::PI_2 * p(d2));
808  break;
809  default:
810  Assert(false, ExcNotImplemented());
811  }
812  }
813  }
814 
815 
816  template <int dim>
817  void
818  CosineGradFunction<dim>::vector_value_list(
819  const std::vector<Point<dim>> &points,
820  std::vector<Vector<double>> & values) const
821  {
822  Assert(values.size() == points.size(),
823  ExcDimensionMismatch(values.size(), points.size()));
824 
825  for (unsigned int i = 0; i < points.size(); ++i)
826  {
827  const Point<dim> &p = points[i];
828  switch (dim)
829  {
830  case 1:
831  values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0));
832  break;
833  case 2:
834  values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
835  std::cos(numbers::PI_2 * p(1));
836  values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
837  std::sin(numbers::PI_2 * p(1));
838  break;
839  case 3:
840  values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
841  std::cos(numbers::PI_2 * p(1)) *
842  std::cos(numbers::PI_2 * p(2));
843  values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
844  std::sin(numbers::PI_2 * p(1)) *
845  std::cos(numbers::PI_2 * p(2));
846  values[i](2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
847  std::cos(numbers::PI_2 * p(1)) *
848  std::sin(numbers::PI_2 * p(2));
849  break;
850  default:
851  Assert(false, ExcNotImplemented());
852  }
853  }
854  }
855 
856 
857  template <int dim>
858  double
859  CosineGradFunction<dim>::laplacian(const Point<dim> & p,
860  const unsigned int d) const
861  {
862  return -numbers::PI_2 * numbers::PI_2 * value(p, d);
863  }
864 
865 
866  template <int dim>
867  Tensor<1, dim>
868  CosineGradFunction<dim>::gradient(const Point<dim> & p,
869  const unsigned int d) const
870  {
871  AssertIndexRange(d, dim);
872  const unsigned int d1 = (d + 1) % dim;
873  const unsigned int d2 = (d + 2) % dim;
874  const double pi2 = numbers::PI_2 * numbers::PI_2;
875 
876  Tensor<1, dim> result;
877  switch (dim)
878  {
879  case 1:
880  result[0] = -pi2 * std::cos(numbers::PI_2 * p(0));
881  break;
882  case 2:
883  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
884  std::cos(numbers::PI_2 * p(d1));
885  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
886  std::sin(numbers::PI_2 * p(d1));
887  break;
888  case 3:
889  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
890  std::cos(numbers::PI_2 * p(d1)) *
891  std::cos(numbers::PI_2 * p(d2));
892  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
893  std::sin(numbers::PI_2 * p(d1)) *
894  std::cos(numbers::PI_2 * p(d2));
895  result[d2] = pi2 * std::sin(numbers::PI_2 * p(d)) *
896  std::cos(numbers::PI_2 * p(d1)) *
897  std::sin(numbers::PI_2 * p(d2));
898  break;
899  default:
900  Assert(false, ExcNotImplemented());
901  }
902  return result;
903  }
904 
905 
906  template <int dim>
907  void
908  CosineGradFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
909  std::vector<Tensor<1, dim>> &gradients,
910  const unsigned int d) const
911  {
912  AssertIndexRange(d, dim);
913  const unsigned int d1 = (d + 1) % dim;
914  const unsigned int d2 = (d + 2) % dim;
915  const double pi2 = numbers::PI_2 * numbers::PI_2;
916 
917  Assert(gradients.size() == points.size(),
918  ExcDimensionMismatch(gradients.size(), points.size()));
919  for (unsigned int i = 0; i < points.size(); ++i)
920  {
921  const Point<dim> &p = points[i];
922  Tensor<1, dim> & result = gradients[i];
923 
924  switch (dim)
925  {
926  case 1:
927  result[0] = -pi2 * std::cos(numbers::PI_2 * p(0));
928  break;
929  case 2:
930  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
931  std::cos(numbers::PI_2 * p(d1));
932  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
933  std::sin(numbers::PI_2 * p(d1));
934  break;
935  case 3:
936  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
937  std::cos(numbers::PI_2 * p(d1)) *
938  std::cos(numbers::PI_2 * p(d2));
939  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
940  std::sin(numbers::PI_2 * p(d1)) *
941  std::cos(numbers::PI_2 * p(d2));
942  result[d2] = pi2 * std::sin(numbers::PI_2 * p(d)) *
943  std::cos(numbers::PI_2 * p(d1)) *
944  std::sin(numbers::PI_2 * p(d2));
945  break;
946  default:
947  Assert(false, ExcNotImplemented());
948  }
949  }
950  }
951 
952 
953  template <int dim>
954  void
955  CosineGradFunction<dim>::vector_gradient_list(
956  const std::vector<Point<dim>> & points,
957  std::vector<std::vector<Tensor<1, dim>>> &gradients) const
958  {
959  AssertVectorVectorDimension(gradients, points.size(), dim);
960  const double pi2 = numbers::PI_2 * numbers::PI_2;
961 
962  for (unsigned int i = 0; i < points.size(); ++i)
963  {
964  const Point<dim> &p = points[i];
965  switch (dim)
966  {
967  case 1:
968  gradients[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p(0));
969  break;
970  case 2:
971  {
972  const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) *
973  std::cos(numbers::PI_2 * p(1));
974  const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) *
975  std::sin(numbers::PI_2 * p(1));
976  gradients[i][0][0] = coco;
977  gradients[i][1][1] = coco;
978  gradients[i][0][1] = sisi;
979  gradients[i][1][0] = sisi;
980  }
981  break;
982  case 3:
983  {
984  const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) *
985  std::cos(numbers::PI_2 * p(1)) *
986  std::cos(numbers::PI_2 * p(2));
987  const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) *
988  std::sin(numbers::PI_2 * p(1)) *
989  std::cos(numbers::PI_2 * p(2));
990  const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) *
991  std::cos(numbers::PI_2 * p(1)) *
992  std::sin(numbers::PI_2 * p(2));
993  const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) *
994  std::sin(numbers::PI_2 * p(1)) *
995  std::sin(numbers::PI_2 * p(2));
996 
997  gradients[i][0][0] = cococo;
998  gradients[i][1][1] = cococo;
999  gradients[i][2][2] = cococo;
1000  gradients[i][0][1] = sisico;
1001  gradients[i][1][0] = sisico;
1002  gradients[i][0][2] = sicosi;
1003  gradients[i][2][0] = sicosi;
1004  gradients[i][1][2] = cosisi;
1005  gradients[i][2][1] = cosisi;
1006  }
1007  break;
1008  default:
1009  Assert(false, ExcNotImplemented());
1010  }
1011  }
1012  }
1013 
1014 
1016 
1017  template <int dim>
1018  double
1019  ExpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1020  {
1021  switch (dim)
1022  {
1023  case 1:
1024  return std::exp(p(0));
1025  case 2:
1026  return std::exp(p(0)) * std::exp(p(1));
1027  case 3:
1028  return std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1029  default:
1030  Assert(false, ExcNotImplemented());
1031  }
1032  return 0.;
1033  }
1034 
1035  template <int dim>
1036  void
1037  ExpFunction<dim>::value_list(const std::vector<Point<dim>> &points,
1038  std::vector<double> & values,
1039  const unsigned int) const
1040  {
1041  Assert(values.size() == points.size(),
1042  ExcDimensionMismatch(values.size(), points.size()));
1043 
1044  for (unsigned int i = 0; i < points.size(); ++i)
1045  {
1046  const Point<dim> &p = points[i];
1047  switch (dim)
1048  {
1049  case 1:
1050  values[i] = std::exp(p(0));
1051  break;
1052  case 2:
1053  values[i] = std::exp(p(0)) * std::exp(p(1));
1054  break;
1055  case 3:
1056  values[i] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1057  break;
1058  default:
1059  Assert(false, ExcNotImplemented());
1060  }
1061  }
1062  }
1063 
1064  template <int dim>
1065  double
1066  ExpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1067  {
1068  switch (dim)
1069  {
1070  case 1:
1071  return std::exp(p(0));
1072  case 2:
1073  return 2 * std::exp(p(0)) * std::exp(p(1));
1074  case 3:
1075  return 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1076  default:
1077  Assert(false, ExcNotImplemented());
1078  }
1079  return 0.;
1080  }
1081 
1082  template <int dim>
1083  void
1084  ExpFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
1085  std::vector<double> & values,
1086  const unsigned int) const
1087  {
1088  Assert(values.size() == points.size(),
1089  ExcDimensionMismatch(values.size(), points.size()));
1090 
1091  for (unsigned int i = 0; i < points.size(); ++i)
1092  {
1093  const Point<dim> &p = points[i];
1094  switch (dim)
1095  {
1096  case 1:
1097  values[i] = std::exp(p(0));
1098  break;
1099  case 2:
1100  values[i] = 2 * std::exp(p(0)) * std::exp(p(1));
1101  break;
1102  case 3:
1103  values[i] = 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1104  break;
1105  default:
1106  Assert(false, ExcNotImplemented());
1107  }
1108  }
1109  }
1110 
1111  template <int dim>
1112  Tensor<1, dim>
1113  ExpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1114  {
1115  Tensor<1, dim> result;
1116  switch (dim)
1117  {
1118  case 1:
1119  result[0] = std::exp(p(0));
1120  break;
1121  case 2:
1122  result[0] = std::exp(p(0)) * std::exp(p(1));
1123  result[1] = result[0];
1124  break;
1125  case 3:
1126  result[0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1127  result[1] = result[0];
1128  result[2] = result[0];
1129  break;
1130  default:
1131  Assert(false, ExcNotImplemented());
1132  }
1133  return result;
1134  }
1135 
1136  template <int dim>
1137  void
1138  ExpFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
1139  std::vector<Tensor<1, dim>> & gradients,
1140  const unsigned int) const
1141  {
1142  Assert(gradients.size() == points.size(),
1143  ExcDimensionMismatch(gradients.size(), points.size()));
1144 
1145  for (unsigned int i = 0; i < points.size(); ++i)
1146  {
1147  const Point<dim> &p = points[i];
1148  switch (dim)
1149  {
1150  case 1:
1151  gradients[i][0] = std::exp(p(0));
1152  break;
1153  case 2:
1154  gradients[i][0] = std::exp(p(0)) * std::exp(p(1));
1155  gradients[i][1] = gradients[i][0];
1156  break;
1157  case 3:
1158  gradients[i][0] =
1159  std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1160  gradients[i][1] = gradients[i][0];
1161  gradients[i][2] = gradients[i][0];
1162  break;
1163  default:
1164  Assert(false, ExcNotImplemented());
1165  }
1166  }
1167  }
1168 
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
1215  LSingularityFunction::vector_value_list(
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
1250  LSingularityFunction::laplacian_list(const std::vector<Point<2>> &points,
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
1308  LSingularityFunction::vector_gradient_list(
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 
1336  LSingularityGradFunction::LSingularityGradFunction()
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
1360  LSingularityGradFunction::value_list(const std::vector<Point<2>> &points,
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
1385  LSingularityGradFunction::vector_value_list(
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
1412  LSingularityGradFunction::laplacian(const Point<2> &,
1413  const unsigned int) const
1414  {
1415  return 0.;
1416  }
1417 
1418 
1419  void
1420  LSingularityGradFunction::laplacian_list(const std::vector<Point<2>> &points,
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>
1434  LSingularityGradFunction::gradient(const Point<2> & /*p*/,
1435  const unsigned int /*component*/) const
1436  {
1437  Assert(false, ExcNotImplemented());
1438  return {};
1439  }
1440 
1441 
1442  void
1443  LSingularityGradFunction::gradient_list(
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
1453  LSingularityGradFunction::vector_gradient_list(
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
1464  SlitSingularityFunction<dim>::value(const Point<dim> &p,
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
1479  SlitSingularityFunction<dim>::value_list(
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
1502  SlitSingularityFunction<dim>::vector_value_list(
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
1527  SlitSingularityFunction<dim>::laplacian(const Point<dim> &,
1528  const unsigned int) const
1529  {
1530  return 0.;
1531  }
1532 
1533 
1534  template <int dim>
1535  void
1536  SlitSingularityFunction<dim>::laplacian_list(
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>
1550  Tensor<1, dim>
1551  SlitSingularityFunction<dim>::gradient(const Point<dim> &p,
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
1572  SlitSingularityFunction<dim>::gradient_list(
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
1601  SlitSingularityFunction<dim>::vector_gradient_list(
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
1634  SlitHyperSingularityFunction::value(const Point<2> &p,
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
1648  SlitHyperSingularityFunction::value_list(const std::vector<Point<2>> &points,
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
1669  SlitHyperSingularityFunction::vector_value_list(
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
1693  SlitHyperSingularityFunction::laplacian(const Point<2> &,
1694  const unsigned int) const
1695  {
1696  return 0.;
1697  }
1698 
1699 
1700  void
1701  SlitHyperSingularityFunction::laplacian_list(
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>
1715  SlitHyperSingularityFunction::gradient(const Point<2> &p,
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
1736  SlitHyperSingularityFunction::gradient_list(
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
1763  SlitHyperSingularityFunction::vector_gradient_list(
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>
1793  JumpFunction<dim>::JumpFunction(const Point<dim> &direction,
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
1827  JumpFunction<dim>::value_list(const std::vector<Point<dim>> &p,
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
1854  JumpFunction<dim>::laplacian_list(const std::vector<Point<dim>> &p,
1855  std::vector<double> & values,
1856  const unsigned int) const
1857  {
1858  Assert(values.size() == p.size(),
1859  ExcDimensionMismatch(values.size(), p.size()));
1860 
1861  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>
1874  Tensor<1, 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
1889  JumpFunction<dim>::gradient_list(const std::vector<Point<dim>> &p,
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
1909  JumpFunction<dim>::memory_consumption() const
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>
1922  FourierCosineFunction<dim>::FourierCosineFunction(
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
1932  FourierCosineFunction<dim>::value(const Point<dim> & p,
1933  const unsigned int component) const
1934  {
1935  (void)component;
1936  Assert(component == 0, ExcIndexRange(component, 0, 1));
1937  return std::cos(fourier_coefficients * p);
1938  }
1939 
1940 
1941 
1942  template <int dim>
1943  Tensor<1, dim>
1944  FourierCosineFunction<dim>::gradient(const Point<dim> & p,
1945  const unsigned int component) const
1946  {
1947  (void)component;
1948  Assert(component == 0, ExcIndexRange(component, 0, 1));
1949  return -fourier_coefficients * std::sin(fourier_coefficients * p);
1950  }
1951 
1952 
1953 
1954  template <int dim>
1955  double
1956  FourierCosineFunction<dim>::laplacian(const Point<dim> & p,
1957  const unsigned int component) const
1958  {
1959  (void)component;
1960  Assert(component == 0, ExcIndexRange(component, 0, 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>
1972  FourierSineFunction<dim>::FourierSineFunction(
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
1982  FourierSineFunction<dim>::value(const Point<dim> & p,
1983  const unsigned int component) const
1984  {
1985  (void)component;
1986  Assert(component == 0, ExcIndexRange(component, 0, 1));
1987  return std::sin(fourier_coefficients * p);
1988  }
1989 
1990 
1991 
1992  template <int dim>
1993  Tensor<1, dim>
1994  FourierSineFunction<dim>::gradient(const Point<dim> & p,
1995  const unsigned int component) const
1996  {
1997  (void)component;
1998  Assert(component == 0, ExcIndexRange(component, 0, 1));
1999  return fourier_coefficients * std::cos(fourier_coefficients * p);
2000  }
2001 
2002 
2003 
2004  template <int dim>
2005  double
2006  FourierSineFunction<dim>::laplacian(const Point<dim> & p,
2007  const unsigned int component) const
2008  {
2009  (void)component;
2010  Assert(component == 0, ExcIndexRange(component, 0, 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>
2022  FourierSineSum<dim>::FourierSineSum(
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(),
2031  ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2032  }
2033 
2034 
2035 
2036  template <int dim>
2037  double
2038  FourierSineSum<dim>::value(const Point<dim> & p,
2039  const unsigned int component) const
2040  {
2041  (void)component;
2042  Assert(component == 0, ExcIndexRange(component, 0, 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>
2055  Tensor<1, dim>
2056  FourierSineSum<dim>::gradient(const Point<dim> & p,
2057  const unsigned int component) const
2058  {
2059  (void)component;
2060  Assert(component == 0, ExcIndexRange(component, 0, 1));
2061 
2062  const unsigned int n = weights.size();
2063  Tensor<1, dim> sum;
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
2074  FourierSineSum<dim>::laplacian(const Point<dim> & p,
2075  const unsigned int component) const
2076  {
2077  (void)component;
2078  Assert(component == 0, ExcIndexRange(component, 0, 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>
2096  FourierCosineSum<dim>::FourierCosineSum(
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(),
2105  ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2106  }
2107 
2108 
2109 
2110  template <int dim>
2111  double
2112  FourierCosineSum<dim>::value(const Point<dim> & p,
2113  const unsigned int component) const
2114  {
2115  (void)component;
2116  Assert(component == 0, ExcIndexRange(component, 0, 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>
2129  Tensor<1, dim>
2130  FourierCosineSum<dim>::gradient(const Point<dim> & p,
2131  const unsigned int component) const
2132  {
2133  (void)component;
2134  Assert(component == 0, ExcIndexRange(component, 0, 1));
2135 
2136  const unsigned int n = weights.size();
2137  Tensor<1, dim> sum;
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
2148  FourierCosineSum<dim>::laplacian(const Point<dim> & p,
2149  const unsigned int component) const
2150  {
2151  (void)component;
2152  Assert(component == 0, ExcIndexRange(component, 0, 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>
2170  Monomial<dim>::Monomial(const Tensor<1, dim> &exponents,
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  Assert(component < this->n_components,
2184  ExcIndexRange(component, 0, this->n_components));
2185 
2186  double prod = 1;
2187  for (unsigned int s = 0; s < dim; ++s)
2188  {
2189  if (p[s] < 0)
2190  Assert(std::floor(exponents[s]) == exponents[s],
2191  ExcMessage("Exponentiation of a negative base number with "
2192  "a real exponent can't be performed."));
2193  prod *= std::pow(p[s], exponents[s]);
2194  }
2195  return prod;
2196  }
2197 
2198 
2199 
2200  template <int dim>
2201  void
2202  Monomial<dim>::vector_value(const Point<dim> &p, Vector<double> &values) const
2203  {
2204  Assert(values.size() == this->n_components,
2205  ExcDimensionMismatch(values.size(), this->n_components));
2206 
2207  for (unsigned int i = 0; i < values.size(); ++i)
2208  values(i) = Monomial<dim>::value(p, i);
2209  }
2210 
2211 
2212 
2213  template <int dim>
2214  Tensor<1, dim>
2215  Monomial<dim>::gradient(const Point<dim> & p,
2216  const unsigned int component) const
2217  {
2218  (void)component;
2219  Assert(component == 0, ExcIndexRange(component, 0, 1));
2220 
2221  Tensor<1, dim> r;
2222  for (unsigned int d = 0; d < dim; ++d)
2223  {
2224  double prod = 1;
2225  for (unsigned int s = 0; s < dim; ++s)
2226  {
2227  if ((s == d) && (exponents[s] == 0) && (p[s] == 0))
2228  {
2229  prod = 0;
2230  break;
2231  }
2232  else
2233  {
2234  if (p[s] < 0)
2235  Assert(std::floor(exponents[s]) == exponents[s],
2236  ExcMessage(
2237  "Exponentiation of a negative base number with "
2238  "a real exponent can't be performed."));
2239  prod *=
2240  (s == d ? exponents[s] * std::pow(p[s], exponents[s] - 1) :
2241  std::pow(p[s], exponents[s]));
2242  }
2243  }
2244  r[d] = prod;
2245  }
2246 
2247  return r;
2248  }
2249 
2250 
2251 
2252  template <int dim>
2253  void
2254  Monomial<dim>::value_list(const std::vector<Point<dim>> &points,
2255  std::vector<double> & values,
2256  const unsigned int component) const
2257  {
2258  Assert(values.size() == points.size(),
2259  ExcDimensionMismatch(values.size(), points.size()));
2260 
2261  for (unsigned int i = 0; i < points.size(); ++i)
2262  values[i] = Monomial<dim>::value(points[i], component);
2263  }
2264 
2265 
2266  template <int dim>
2267  Bessel1<dim>::Bessel1(const unsigned int order,
2268  const double wave_number,
2269  const Point<dim> center)
2270  : order(order)
2271  , wave_number(wave_number)
2272  , center(center)
2273  {
2274  Assert(wave_number >= 0., ExcMessage("wave_number must be nonnegative!"));
2275  }
2276 
2277  template <int dim>
2278  double
2279  Bessel1<dim>::value(const Point<dim> &p, const unsigned int) const
2280  {
2281  Assert(dim == 2, ExcNotImplemented());
2282  const double r = p.distance(center);
2283  return std_cxx17::cyl_bessel_j(order, r * wave_number);
2284  }
2285 
2286 
2287  template <int dim>
2288  void
2289  Bessel1<dim>::value_list(const std::vector<Point<dim>> &points,
2290  std::vector<double> & values,
2291  const unsigned int) const
2292  {
2293  Assert(dim == 2, ExcNotImplemented());
2294  AssertDimension(points.size(), values.size());
2295  for (unsigned int k = 0; k < points.size(); ++k)
2296  {
2297  const double r = points[k].distance(center);
2298  values[k] = std_cxx17::cyl_bessel_j(order, r * wave_number);
2299  }
2300  }
2301 
2302 
2303  template <int dim>
2304  Tensor<1, dim>
2305  Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
2306  {
2307  Assert(dim == 2, ExcNotImplemented());
2308  const double r = p.distance(center);
2309  const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
2310  const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
2311 
2312  const double dJn =
2313  (order == 0) ?
2314  (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2315  (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2316  std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2317  Tensor<1, dim> result;
2318  result[0] = wave_number * co * dJn;
2319  result[1] = wave_number * si * dJn;
2320  return result;
2321  }
2322 
2323 
2324 
2325  template <int dim>
2326  void
2327  Bessel1<dim>::gradient_list(const std::vector<Point<dim>> &points,
2328  std::vector<Tensor<1, dim>> & gradients,
2329  const unsigned int) const
2330  {
2331  Assert(dim == 2, ExcNotImplemented());
2332  AssertDimension(points.size(), gradients.size());
2333  for (unsigned int k = 0; k < points.size(); ++k)
2334  {
2335  const Point<dim> &p = points[k];
2336  const double r = p.distance(center);
2337  const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
2338  const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
2339 
2340  const double dJn =
2341  (order == 0) ?
2342  (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2343  (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2344  std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2345  Tensor<1, dim> &result = gradients[k];
2346  result[0] = wave_number * co * dJn;
2347  result[1] = wave_number * si * dJn;
2348  }
2349  }
2350 
2351 
2352 
2353  namespace
2354  {
2355  // interpolate a data value from a table where ix denotes
2356  // the (lower) left endpoint of the interval to interpolate
2357  // in, and p_unit denotes the point in unit coordinates to do so.
2358  double
2359  interpolate(const Table<1, double> &data_values,
2360  const TableIndices<1> & ix,
2361  const Point<1> & xi)
2362  {
2363  return ((1 - xi[0]) * data_values[ix[0]] +
2364  xi[0] * data_values[ix[0] + 1]);
2365  }
2366 
2367  double
2368  interpolate(const Table<2, double> &data_values,
2369  const TableIndices<2> & ix,
2370  const Point<2> & p_unit)
2371  {
2372  return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
2373  p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
2374  (1 - p_unit[1]) +
2375  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
2376  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
2377  p_unit[1]);
2378  }
2379 
2380  double
2381  interpolate(const Table<3, double> &data_values,
2382  const TableIndices<3> & ix,
2383  const Point<3> & p_unit)
2384  {
2385  return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
2386  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
2387  (1 - p_unit[1]) +
2388  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
2389  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
2390  p_unit[1]) *
2391  (1 - p_unit[2]) +
2392  (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
2393  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
2394  (1 - p_unit[1]) +
2395  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
2396  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
2397  p_unit[1]) *
2398  p_unit[2]);
2399  }
2400 
2401 
2402  // Interpolate the gradient of a data value from a table where ix
2403  // denotes the lower left endpoint of the interval to interpolate
2404  // in, p_unit denotes the point in unit coordinates, and dx
2405  // denotes the width of the interval in each dimension.
2406  Tensor<1, 1>
2407  gradient_interpolate(const Table<1, double> &data_values,
2408  const TableIndices<1> & ix,
2409  const Point<1> & p_unit,
2410  const Point<1> & dx)
2411  {
2412  (void)p_unit;
2413  Tensor<1, 1> grad;
2414  grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
2415  return grad;
2416  }
2417 
2418 
2419  Tensor<1, 2>
2420  gradient_interpolate(const Table<2, double> &data_values,
2421  const TableIndices<2> & ix,
2422  const Point<2> & p_unit,
2423  const Point<2> & dx)
2424  {
2425  Tensor<1, 2> grad;
2426  double u00 = data_values[ix[0]][ix[1]],
2427  u01 = data_values[ix[0] + 1][ix[1]],
2428  u10 = data_values[ix[0]][ix[1] + 1],
2429  u11 = data_values[ix[0] + 1][ix[1] + 1];
2430 
2431  grad[0] =
2432  ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
2433  grad[1] =
2434  ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
2435  return grad;
2436  }
2437 
2438 
2439  Tensor<1, 3>
2440  gradient_interpolate(const Table<3, double> &data_values,
2441  const TableIndices<3> & ix,
2442  const Point<3> & p_unit,
2443  const Point<3> & dx)
2444  {
2445  Tensor<1, 3> grad;
2446  double u000 = data_values[ix[0]][ix[1]][ix[2]],
2447  u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
2448  u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
2449  u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
2450  u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
2451  u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
2452  u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
2453  u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
2454 
2455  grad[0] =
2456  ((1 - p_unit[2]) *
2457  ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
2458  p_unit[2] *
2459  ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
2460  dx[0];
2461  grad[1] =
2462  ((1 - p_unit[2]) *
2463  ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
2464  p_unit[2] *
2465  ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
2466  dx[1];
2467  grad[2] =
2468  ((1 - p_unit[1]) *
2469  ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
2470  p_unit[1] *
2471  ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
2472  dx[2];
2473 
2474  return grad;
2475  }
2476  } // namespace
2477 
2478  template <int dim>
2479  InterpolatedTensorProductGridData<dim>::InterpolatedTensorProductGridData(
2480  const std::array<std::vector<double>, dim> &coordinate_values,
2481  const Table<dim, double> & data_values)
2482  : coordinate_values(coordinate_values)
2483  , data_values(data_values)
2484  {
2485  for (unsigned int d = 0; d < dim; ++d)
2486  {
2487  Assert(
2488  coordinate_values[d].size() >= 2,
2489  ExcMessage(
2490  "Coordinate arrays must have at least two coordinate values!"));
2491  for (unsigned int i = 0; i < coordinate_values[d].size() - 1; ++i)
2492  Assert(
2493  coordinate_values[d][i] < coordinate_values[d][i + 1],
2494  ExcMessage(
2495  "Coordinate arrays must be sorted in strictly ascending order."));
2496 
2497  Assert(data_values.size()[d] == coordinate_values[d].size(),
2498  ExcMessage(
2499  "Data and coordinate tables do not have the same size."));
2500  }
2501  }
2502 
2503 
2504 
2505  template <int dim>
2506  TableIndices<dim>
2507  InterpolatedTensorProductGridData<dim>::table_index_of_point(
2508  const Point<dim> &p) const
2509  {
2510  // find out where this data point lies, relative to the given
2511  // points. if we run all the way to the end of the range,
2512  // set the indices so that we will simply query the last of the
2513  // intervals, starting at x.size()-2 and going to x.size()-1.
2514  TableIndices<dim> ix;
2515  for (unsigned int d = 0; d < dim; ++d)
2516  {
2517  // get the index of the first element of the coordinate arrays that is
2518  // larger than p[d]
2519  ix[d] = (std::lower_bound(coordinate_values[d].begin(),
2520  coordinate_values[d].end(),
2521  p[d]) -
2522  coordinate_values[d].begin());
2523 
2524  // the one we want is the index of the coordinate to the left, however,
2525  // so decrease it by one (unless we have a point to the left of all, in
2526  // which case we stay where we are; the formulas below are made in a way
2527  // that allow us to extend the function by a constant value)
2528  //
2529  // to make this work, if we got coordinate_values[d].end(), we actually
2530  // have to consider the last box which has index size()-2
2531  if (ix[d] == coordinate_values[d].size())
2532  ix[d] = coordinate_values[d].size() - 2;
2533  else if (ix[d] > 0)
2534  --ix[d];
2535  }
2536 
2537  return ix;
2538  }
2539 
2540 
2541 
2542  template <int dim>
2543  double
2544  InterpolatedTensorProductGridData<dim>::value(
2545  const Point<dim> & p,
2546  const unsigned int component) const
2547  {
2548  (void)component;
2549  Assert(
2550  component == 0,
2551  ExcMessage(
2552  "This is a scalar function object, the component can only be zero."));
2553 
2554  // find the index in the data table of the cell containing the input point
2555  const TableIndices<dim> ix = table_index_of_point(p);
2556 
2557  // now compute the relative point within the interval/rectangle/box
2558  // defined by the point coordinates found above. truncate below and
2559  // above to accommodate points that may lie outside the range
2560  Point<dim> p_unit;
2561  for (unsigned int d = 0; d < dim; ++d)
2562  p_unit[d] = std::max(std::min((p[d] - coordinate_values[d][ix[d]]) /
2563  (coordinate_values[d][ix[d] + 1] -
2564  coordinate_values[d][ix[d]]),
2565  1.),
2566  0.);
2567 
2568  return interpolate(data_values, ix, p_unit);
2569  }
2570 
2571 
2572 
2573  template <int dim>
2574  Tensor<1, dim>
2575  InterpolatedTensorProductGridData<dim>::gradient(
2576  const Point<dim> & p,
2577  const unsigned int component) const
2578  {
2579  (void)component;
2580  Assert(
2581  component == 0,
2582  ExcMessage(
2583  "This is a scalar function object, the component can only be zero."));
2584 
2585  // find out where this data point lies
2586  const TableIndices<dim> ix = table_index_of_point(p);
2587 
2588  Point<dim> dx;
2589  for (unsigned int d = 0; d < dim; ++d)
2590  dx[d] = coordinate_values[d][ix[d] + 1] - coordinate_values[d][ix[d]];
2591 
2592  Point<dim> p_unit;
2593  for (unsigned int d = 0; d < dim; ++d)
2594  p_unit[d] =
2595  std::max(std::min((p[d] - coordinate_values[d][ix[d]]) / dx[d], 1.),
2596  0.0);
2597 
2598  return gradient_interpolate(data_values, ix, p_unit, dx);
2599  }
2600 
2601 
2602 
2603  template <int dim>
2604  InterpolatedUniformGridData<dim>::InterpolatedUniformGridData(
2605  const std::array<std::pair<double, double>, dim> &interval_endpoints,
2606  const std::array<unsigned int, dim> & n_subintervals,
2607  const Table<dim, double> & data_values)
2608  : interval_endpoints(interval_endpoints)
2609  , n_subintervals(n_subintervals)
2610  , data_values(data_values)
2611  {
2612  for (unsigned int d = 0; d < dim; ++d)
2613  {
2614  Assert(n_subintervals[d] >= 1,
2615  ExcMessage("There needs to be at least one subinterval in each "
2616  "coordinate direction."));
2617  Assert(interval_endpoints[d].first < interval_endpoints[d].second,
2618  ExcMessage("The interval in each coordinate direction needs "
2619  "to have positive size"));
2620  Assert(data_values.size()[d] == n_subintervals[d] + 1,
2621  ExcMessage("The data table does not have the correct size."));
2622  }
2623  }
2624 
2625 
2626  template <int dim>
2627  double
2628  InterpolatedUniformGridData<dim>::value(const Point<dim> & p,
2629  const unsigned int component) const
2630  {
2631  (void)component;
2632  Assert(
2633  component == 0,
2634  ExcMessage(
2635  "This is a scalar function object, the component can only be zero."));
2636 
2637  // find out where this data point lies, relative to the given
2638  // subdivision points
2639  TableIndices<dim> ix;
2640  for (unsigned int d = 0; d < dim; ++d)
2641  {
2642  const double delta_x =
2643  ((interval_endpoints[d].second - interval_endpoints[d].first) /
2644  n_subintervals[d]);
2645  if (p[d] <= interval_endpoints[d].first)
2646  ix[d] = 0;
2647  else if (p[d] >= interval_endpoints[d].second - delta_x)
2648  ix[d] = n_subintervals[d] - 1;
2649  else
2650  ix[d] = static_cast<unsigned int>(
2651  (p[d] - interval_endpoints[d].first) / delta_x);
2652  }
2653 
2654  // now compute the relative point within the interval/rectangle/box
2655  // defined by the point coordinates found above. truncate below and
2656  // above to accommodate points that may lie outside the range
2657  Point<dim> p_unit;
2658  for (unsigned int d = 0; d < dim; ++d)
2659  {
2660  const double delta_x =
2661  ((interval_endpoints[d].second - interval_endpoints[d].first) /
2662  n_subintervals[d]);
2663  p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
2664  ix[d] * delta_x) /
2665  delta_x,
2666  1.),
2667  0.);
2668  }
2669 
2670  return interpolate(data_values, ix, p_unit);
2671  }
2672 
2673  /* ---------------------- Polynomial ----------------------- */
2674 
2675 
2676 
2677  template <int dim>
2678  Polynomial<dim>::Polynomial(const Table<2, double> & exponents,
2679  const std::vector<double> &coefficients)
2680  : Function<dim>(1)
2681  , exponents(exponents)
2682  , coefficients(coefficients)
2683  {
2684  Assert(exponents.n_rows() == coefficients.size(),
2685  ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
2686  Assert(exponents.n_cols() == dim,
2687  ExcDimensionMismatch(exponents.n_cols(), dim));
2688  }
2689 
2690 
2691 
2692  template <int dim>
2693  double
2694  Polynomial<dim>::value(const Point<dim> & p,
2695  const unsigned int component) const
2696  {
2697  (void)component;
2698  Assert(component == 0, ExcIndexRange(component, 0, 1));
2699 
2700  double sum = 0;
2701  for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2702  {
2703  double prod = 1;
2704  for (unsigned int s = 0; s < dim; ++s)
2705  {
2706  if (p[s] < 0)
2707  Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
2708  ExcMessage("Exponentiation of a negative base number with "
2709  "a real exponent can't be performed."));
2710  prod *= std::pow(p[s], exponents[monom][s]);
2711  }
2712  sum += coefficients[monom] * prod;
2713  }
2714  return sum;
2715  }
2716 
2717 
2718 
2719  template <int dim>
2720  void
2721  Polynomial<dim>::value_list(const std::vector<Point<dim>> &points,
2722  std::vector<double> & values,
2723  const unsigned int component) const
2724  {
2725  Assert(values.size() == points.size(),
2726  ExcDimensionMismatch(values.size(), points.size()));
2727 
2728  for (unsigned int i = 0; i < points.size(); ++i)
2729  values[i] = Polynomial<dim>::value(points[i], component);
2730  }
2731 
2732 
2733 
2734  template <int dim>
2735  Tensor<1, dim>
2736  Polynomial<dim>::gradient(const Point<dim> & p,
2737  const unsigned int component) const
2738  {
2739  (void)component;
2740  Assert(component == 0, ExcIndexRange(component, 0, 1));
2741 
2742  Tensor<1, dim> r;
2743 
2744  for (unsigned int d = 0; d < dim; ++d)
2745  {
2746  double sum = 0;
2747 
2748  for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2749  {
2750  double prod = 1;
2751  for (unsigned int s = 0; s < dim; ++s)
2752  {
2753  if ((s == d) && (exponents[monom][s] == 0) && (p[s] == 0))
2754  {
2755  prod = 0;
2756  break;
2757  }
2758  else
2759  {
2760  if (p[s] < 0)
2761  Assert(std::floor(exponents[monom][s]) ==
2762  exponents[monom][s],
2763  ExcMessage(
2764  "Exponentiation of a negative base number with "
2765  "a real exponent can't be performed."));
2766  prod *=
2767  (s == d ? exponents[monom][s] *
2768  std::pow(p[s], exponents[monom][s] - 1) :
2769  std::pow(p[s], exponents[monom][s]));
2770  }
2771  }
2772  sum += coefficients[monom] * prod;
2773  }
2774  r[d] = sum;
2775  }
2776  return r;
2777  }
2778 
2779 
2780  // explicit instantiations
2781  template class SquareFunction<1>;
2782  template class SquareFunction<2>;
2783  template class SquareFunction<3>;
2784  template class Q1WedgeFunction<1>;
2785  template class Q1WedgeFunction<2>;
2786  template class Q1WedgeFunction<3>;
2787  template class PillowFunction<1>;
2788  template class PillowFunction<2>;
2789  template class PillowFunction<3>;
2790  template class CosineFunction<1>;
2791  template class CosineFunction<2>;
2792  template class CosineFunction<3>;
2793  template class CosineGradFunction<1>;
2794  template class CosineGradFunction<2>;
2795  template class CosineGradFunction<3>;
2796  template class ExpFunction<1>;
2797  template class ExpFunction<2>;
2798  template class ExpFunction<3>;
2799  template class JumpFunction<1>;
2800  template class JumpFunction<2>;
2801  template class JumpFunction<3>;
2802  template class FourierCosineFunction<1>;
2803  template class FourierCosineFunction<2>;
2804  template class FourierCosineFunction<3>;
2805  template class FourierSineFunction<1>;
2806  template class FourierSineFunction<2>;
2807  template class FourierSineFunction<3>;
2808  template class FourierCosineSum<1>;
2809  template class FourierCosineSum<2>;
2810  template class FourierCosineSum<3>;
2811  template class FourierSineSum<1>;
2812  template class FourierSineSum<2>;
2813  template class FourierSineSum<3>;
2814  template class SlitSingularityFunction<2>;
2815  template class SlitSingularityFunction<3>;
2816  template class Monomial<1>;
2817  template class Monomial<2>;
2818  template class Monomial<3>;
2819  template class Bessel1<1>;
2820  template class Bessel1<2>;
2821  template class Bessel1<3>;
2822  template class InterpolatedTensorProductGridData<1>;
2823  template class InterpolatedTensorProductGridData<2>;
2824  template class InterpolatedTensorProductGridData<3>;
2825  template class InterpolatedUniformGridData<1>;
2826  template class InterpolatedUniformGridData<2>;
2827  template class InterpolatedUniformGridData<3>;
2828  template class Polynomial<1>;
2829  template class Polynomial<2>;
2830  template class Polynomial<3>;
2831 } // namespace Functions
2832 
2833 DEAL_II_NAMESPACE_CLOSE
Functions::FourierCosineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2130
Functions::CosineGradFunction::value
virtual double value(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:722
Functions::CosineFunction::hessian
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:598
Functions::PillowFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:261
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Functions::InterpolatedTensorProductGridData::table_index_of_point
TableIndices< dim > table_index_of_point(const Point< dim > &p) const
Definition: function_lib.cc:2507
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::PillowFunction::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:280
Functions::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::Monomial::Monomial
Monomial(const Tensor< 1, dim > &exponents, const unsigned int n_components=1)
Definition: function_lib.cc:2170
Functions::Polynomial::Polynomial
Polynomial(const Table< 2, double > &exponents, const std::vector< double > &coefficients)
Definition: function_lib.cc:2678
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
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
Functions::PillowFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:364
Functions::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::Monomial::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:2202
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1595
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::FourierCosineSum::FourierCosineSum
FourierCosineSum(const std::vector< Point< dim >> &fourier_coefficients, const std::vector< double > &weights)
Definition: function_lib.cc:2096
Functions::FourierCosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1956
Functions::SquareFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:34
Functions::FourierSineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1982
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Functions::Monomial::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2180
Functions::InterpolatedTensorProductGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2544
Functions::FourierSineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1994
Functions::CosineFunction::CosineFunction
CosineFunction(const unsigned int n_components=1)
Definition: function_lib.cc:425
Functions::ExpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1066
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::FourierCosineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1932
Functions::Monomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2215
Functions::FourierCosineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2148
Functions::InterpolatedTensorProductGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2575
Functions::Q1WedgeFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:176
Functions::Bessel1::value
virtual double value(const Point< dim > &points, const unsigned int component=0) const override
Definition: function_lib.cc:2279
Functions::SlitSingularityFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1551
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
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::FourierCosineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2112
Functions::JumpFunction::memory_consumption
std::size_t memory_consumption() const
Definition: function_lib.cc:1909
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:2254
Functions::CosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:486
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::JumpFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1875
Functions::SquareFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:42
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
Functions::FourierSineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2074
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Functions::SlitSingularityFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1527
Functions::InterpolatedUniformGridData::data_values
const Table< dim, double > data_values
Definition: function_lib.h:1579
Functions::Bessel1::Bessel1
Bessel1(const unsigned int order, const double wave_number, const Point< dim > center=Point< dim >())
Definition: function_lib.cc:2267
Functions::CosineGradFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:748
Functions::SlitSingularityFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1464
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::FourierCosineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition: function_lib.h:906
numbers::PI_2
static constexpr double PI_2
Definition: numbers.h:151
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
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::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:868
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::SquareFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:69
Functions::Polynomial::exponents
const Table< 2, double > exponents
Definition: function_lib.h:1638
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::Polynomial::coefficients
const std::vector< double > coefficients
Definition: function_lib.h:1643
Functions::Q1WedgeFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:129
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::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::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::JumpFunction::direction
const Point< dim > direction
Definition: function_lib.h:669
TableBase::size
size_type size(const unsigned int i) const
Functions::FourierSineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition: function_lib.h:851
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:2289
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::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::JumpFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1817
Functions::PillowFunction::PillowFunction
PillowFunction(const double offset=0.)
Definition: function_lib.cc:254
numbers::PI
static constexpr double PI
Definition: numbers.h:146
Functions::InterpolatedTensorProductGridData::coordinate_values
const std::array< std::vector< double >, dim > coordinate_values
Definition: function_lib.h:1484
Functions::InterpolatedTensorProductGridData::InterpolatedTensorProductGridData
InterpolatedTensorProductGridData(const std::array< std::vector< double >, dim > &coordinate_values, const Table< dim, double > &data_values)
Definition: function_lib.cc:2479
Functions::InterpolatedUniformGridData::n_subintervals
const std::array< unsigned int, dim > n_subintervals
Definition: function_lib.h:1574
Point::square
numbers::NumberTraits< Number >::real_type square() const
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:2604
Functions::FourierCosineFunction::FourierCosineFunction
FourierCosineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition: function_lib.cc:1922
Vector::size
size_type size() const
Functions::PillowFunction::gradient_list
virtual void gradient_list(const std::vector< Point< dim >> &points, std::vector< Tensor< 1, dim >> &gradients, const unsigned int component=0) const override
Definition: function_lib.cc:389
Functions::ExpFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1019
Functions::PillowFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:313
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Functions::FourierSineSum::FourierSineSum
FourierSineSum(const std::vector< Point< dim >> &fourier_coefficients, const std::vector< double > &weights)
Definition: function_lib.cc:2022
Functions::CosineFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:433
Functions::FourierSineFunction::FourierSineFunction
FourierSineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition: function_lib.cc:1972
Table
Definition: table.h:37
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::InterpolatedUniformGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2628
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::Polynomial::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2694
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
Functions::Polynomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2736
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
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::CosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:523
Functions::InterpolatedUniformGridData::interval_endpoints
const std::array< std::pair< double, double >, dim > interval_endpoints
Definition: function_lib.h:1569
StandardExceptions::ExcZero
static ::ExceptionBase & ExcZero()
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::CosineGradFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:859
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::FourierSineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2038
Functions::JumpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1844
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::FourierSineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2006
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
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
Functions::Bessel1::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2305
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
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:2721
Functions::JumpFunction::JumpFunction
JumpFunction(const Point< dim > &direction, const double steepness)
Definition: function_lib.cc:1793
Generated by   doxygen 1.8.14