Reference documentation for deal.II version 9.1.0-pre
function_lib.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2017 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  if (true)
610  {
611  const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) *
612  std::cos(numbers::PI_2 * p(1));
613  const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) *
614  std::sin(numbers::PI_2 * p(1));
615  result[0][0] = coco;
616  result[1][1] = coco;
617  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
618  result[0][1] = sisi;
619  }
620  break;
621  case 3:
622  if (true)
623  {
624  const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) *
625  std::cos(numbers::PI_2 * p(1)) *
626  std::cos(numbers::PI_2 * p(2));
627  const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) *
628  std::sin(numbers::PI_2 * p(1)) *
629  std::cos(numbers::PI_2 * p(2));
630  const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) *
631  std::cos(numbers::PI_2 * p(1)) *
632  std::sin(numbers::PI_2 * p(2));
633  const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) *
634  std::sin(numbers::PI_2 * p(1)) *
635  std::sin(numbers::PI_2 * p(2));
636 
637  result[0][0] = cococo;
638  result[1][1] = cococo;
639  result[2][2] = cococo;
640  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
641  result[0][1] = sisico;
642  result[0][2] = sicosi;
643  result[1][2] = cosisi;
644  }
645  break;
646  default:
647  Assert(false, ExcNotImplemented());
648  }
649  return result;
650  }
651 
652  template <int dim>
653  void
654  CosineFunction<dim>::hessian_list(
655  const std::vector<Point<dim>> & points,
656  std::vector<SymmetricTensor<2, dim>> &hessians,
657  const unsigned int) const
658  {
659  Assert(hessians.size() == points.size(),
660  ExcDimensionMismatch(hessians.size(), points.size()));
661 
662  const double pi2 = numbers::PI_2 * numbers::PI_2;
663 
664  for (unsigned int i = 0; i < points.size(); ++i)
665  {
666  const Point<dim> &p = points[i];
667  switch (dim)
668  {
669  case 1:
670  hessians[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p(0));
671  break;
672  case 2:
673  if (true)
674  {
675  const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) *
676  std::cos(numbers::PI_2 * p(1));
677  const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) *
678  std::sin(numbers::PI_2 * p(1));
679  hessians[i][0][0] = coco;
680  hessians[i][1][1] = coco;
681  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
682  hessians[i][0][1] = sisi;
683  }
684  break;
685  case 3:
686  if (true)
687  {
688  const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) *
689  std::cos(numbers::PI_2 * p(1)) *
690  std::cos(numbers::PI_2 * p(2));
691  const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) *
692  std::sin(numbers::PI_2 * p(1)) *
693  std::cos(numbers::PI_2 * p(2));
694  const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) *
695  std::cos(numbers::PI_2 * p(1)) *
696  std::sin(numbers::PI_2 * p(2));
697  const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) *
698  std::sin(numbers::PI_2 * p(1)) *
699  std::sin(numbers::PI_2 * p(2));
700 
701  hessians[i][0][0] = cococo;
702  hessians[i][1][1] = cococo;
703  hessians[i][2][2] = cococo;
704  // for SymmetricTensor we assign [ij] and [ji] simultaneously:
705  hessians[i][0][1] = sisico;
706  hessians[i][0][2] = sicosi;
707  hessians[i][1][2] = cosisi;
708  }
709  break;
710  default:
711  Assert(false, ExcNotImplemented());
712  }
713  }
714  }
715 
717 
718  template <int dim>
719  CosineGradFunction<dim>::CosineGradFunction()
720  : Function<dim>(dim)
721  {}
722 
723 
724  template <int dim>
725  double
726  CosineGradFunction<dim>::value(const Point<dim> & p,
727  const unsigned int d) const
728  {
729  AssertIndexRange(d, dim);
730  const unsigned int d1 = (d + 1) % dim;
731  const unsigned int d2 = (d + 2) % dim;
732  switch (dim)
733  {
734  case 1:
735  return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(0)));
736  case 2:
737  return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
738  std::cos(numbers::PI_2 * p(d1)));
739  case 3:
740  return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
741  std::cos(numbers::PI_2 * p(d1)) *
742  std::cos(numbers::PI_2 * p(d2)));
743  default:
744  Assert(false, ExcNotImplemented());
745  }
746  return 0.;
747  }
748 
749 
750  template <int dim>
751  void
752  CosineGradFunction<dim>::vector_value(const Point<dim> &p,
753  Vector<double> & result) const
754  {
755  AssertDimension(result.size(), dim);
756  switch (dim)
757  {
758  case 1:
759  result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0));
760  break;
761  case 2:
762  result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
763  std::cos(numbers::PI_2 * p(1));
764  result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
765  std::sin(numbers::PI_2 * p(1));
766  break;
767  case 3:
768  result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
769  std::cos(numbers::PI_2 * p(1)) *
770  std::cos(numbers::PI_2 * p(2));
771  result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
772  std::sin(numbers::PI_2 * p(1)) *
773  std::cos(numbers::PI_2 * p(2));
774  result(2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
775  std::cos(numbers::PI_2 * p(1)) *
776  std::sin(numbers::PI_2 * p(2));
777  break;
778  default:
779  Assert(false, ExcNotImplemented());
780  }
781  }
782 
783 
784  template <int dim>
785  void
786  CosineGradFunction<dim>::value_list(const std::vector<Point<dim>> &points,
787  std::vector<double> & values,
788  const unsigned int d) const
789  {
790  Assert(values.size() == points.size(),
791  ExcDimensionMismatch(values.size(), points.size()));
792  AssertIndexRange(d, dim);
793  const unsigned int d1 = (d + 1) % dim;
794  const unsigned int d2 = (d + 2) % dim;
795 
796  for (unsigned int i = 0; i < points.size(); ++i)
797  {
798  const Point<dim> &p = points[i];
799  switch (dim)
800  {
801  case 1:
802  values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d));
803  break;
804  case 2:
805  values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
806  std::cos(numbers::PI_2 * p(d1));
807  break;
808  case 3:
809  values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) *
810  std::cos(numbers::PI_2 * p(d1)) *
811  std::cos(numbers::PI_2 * p(d2));
812  break;
813  default:
814  Assert(false, ExcNotImplemented());
815  }
816  }
817  }
818 
819 
820  template <int dim>
821  void
822  CosineGradFunction<dim>::vector_value_list(
823  const std::vector<Point<dim>> &points,
824  std::vector<Vector<double>> & values) const
825  {
826  Assert(values.size() == points.size(),
827  ExcDimensionMismatch(values.size(), points.size()));
828 
829  for (unsigned int i = 0; i < points.size(); ++i)
830  {
831  const Point<dim> &p = points[i];
832  switch (dim)
833  {
834  case 1:
835  values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0));
836  break;
837  case 2:
838  values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
839  std::cos(numbers::PI_2 * p(1));
840  values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
841  std::sin(numbers::PI_2 * p(1));
842  break;
843  case 3:
844  values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) *
845  std::cos(numbers::PI_2 * p(1)) *
846  std::cos(numbers::PI_2 * p(2));
847  values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
848  std::sin(numbers::PI_2 * p(1)) *
849  std::cos(numbers::PI_2 * p(2));
850  values[i](2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) *
851  std::cos(numbers::PI_2 * p(1)) *
852  std::sin(numbers::PI_2 * p(2));
853  break;
854  default:
855  Assert(false, ExcNotImplemented());
856  }
857  }
858  }
859 
860 
861  template <int dim>
862  double
863  CosineGradFunction<dim>::laplacian(const Point<dim> & p,
864  const unsigned int d) const
865  {
866  return -numbers::PI_2 * numbers::PI_2 * value(p, d);
867  }
868 
869 
870  template <int dim>
871  Tensor<1, dim>
872  CosineGradFunction<dim>::gradient(const Point<dim> & p,
873  const unsigned int d) const
874  {
875  AssertIndexRange(d, dim);
876  const unsigned int d1 = (d + 1) % dim;
877  const unsigned int d2 = (d + 2) % dim;
878  const double pi2 = numbers::PI_2 * numbers::PI_2;
879 
880  Tensor<1, dim> result;
881  switch (dim)
882  {
883  case 1:
884  result[0] = -pi2 * std::cos(numbers::PI_2 * p(0));
885  break;
886  case 2:
887  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
888  std::cos(numbers::PI_2 * p(d1));
889  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
890  std::sin(numbers::PI_2 * p(d1));
891  break;
892  case 3:
893  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
894  std::cos(numbers::PI_2 * p(d1)) *
895  std::cos(numbers::PI_2 * p(d2));
896  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
897  std::sin(numbers::PI_2 * p(d1)) *
898  std::cos(numbers::PI_2 * p(d2));
899  result[d2] = pi2 * std::sin(numbers::PI_2 * p(d)) *
900  std::cos(numbers::PI_2 * p(d1)) *
901  std::sin(numbers::PI_2 * p(d2));
902  break;
903  default:
904  Assert(false, ExcNotImplemented());
905  }
906  return result;
907  }
908 
909 
910  template <int dim>
911  void
912  CosineGradFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
913  std::vector<Tensor<1, dim>> &gradients,
914  const unsigned int d) const
915  {
916  AssertIndexRange(d, dim);
917  const unsigned int d1 = (d + 1) % dim;
918  const unsigned int d2 = (d + 2) % dim;
919  const double pi2 = numbers::PI_2 * numbers::PI_2;
920 
921  Assert(gradients.size() == points.size(),
922  ExcDimensionMismatch(gradients.size(), points.size()));
923  for (unsigned int i = 0; i < points.size(); ++i)
924  {
925  const Point<dim> &p = points[i];
926  Tensor<1, dim> & result = gradients[i];
927 
928  switch (dim)
929  {
930  case 1:
931  result[0] = -pi2 * std::cos(numbers::PI_2 * p(0));
932  break;
933  case 2:
934  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
935  std::cos(numbers::PI_2 * p(d1));
936  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
937  std::sin(numbers::PI_2 * p(d1));
938  break;
939  case 3:
940  result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) *
941  std::cos(numbers::PI_2 * p(d1)) *
942  std::cos(numbers::PI_2 * p(d2));
943  result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) *
944  std::sin(numbers::PI_2 * p(d1)) *
945  std::cos(numbers::PI_2 * p(d2));
946  result[d2] = pi2 * std::sin(numbers::PI_2 * p(d)) *
947  std::cos(numbers::PI_2 * p(d1)) *
948  std::sin(numbers::PI_2 * p(d2));
949  break;
950  default:
951  Assert(false, ExcNotImplemented());
952  }
953  }
954  }
955 
956 
957  template <int dim>
958  void
959  CosineGradFunction<dim>::vector_gradient_list(
960  const std::vector<Point<dim>> & points,
961  std::vector<std::vector<Tensor<1, dim>>> &gradients) const
962  {
963  AssertVectorVectorDimension(gradients, points.size(), dim);
964  const double pi2 = numbers::PI_2 * numbers::PI_2;
965 
966  for (unsigned int i = 0; i < points.size(); ++i)
967  {
968  const Point<dim> &p = points[i];
969  switch (dim)
970  {
971  case 1:
972  gradients[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p(0));
973  break;
974  case 2:
975  if (true)
976  {
977  const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) *
978  std::cos(numbers::PI_2 * p(1));
979  const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) *
980  std::sin(numbers::PI_2 * p(1));
981  gradients[i][0][0] = coco;
982  gradients[i][1][1] = coco;
983  gradients[i][0][1] = sisi;
984  gradients[i][1][0] = sisi;
985  }
986  break;
987  case 3:
988  if (true)
989  {
990  const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) *
991  std::cos(numbers::PI_2 * p(1)) *
992  std::cos(numbers::PI_2 * p(2));
993  const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) *
994  std::sin(numbers::PI_2 * p(1)) *
995  std::cos(numbers::PI_2 * p(2));
996  const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) *
997  std::cos(numbers::PI_2 * p(1)) *
998  std::sin(numbers::PI_2 * p(2));
999  const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) *
1000  std::sin(numbers::PI_2 * p(1)) *
1001  std::sin(numbers::PI_2 * p(2));
1002 
1003  gradients[i][0][0] = cococo;
1004  gradients[i][1][1] = cococo;
1005  gradients[i][2][2] = cococo;
1006  gradients[i][0][1] = sisico;
1007  gradients[i][1][0] = sisico;
1008  gradients[i][0][2] = sicosi;
1009  gradients[i][2][0] = sicosi;
1010  gradients[i][1][2] = cosisi;
1011  gradients[i][2][1] = cosisi;
1012  }
1013  break;
1014  default:
1015  Assert(false, ExcNotImplemented());
1016  }
1017  }
1018  }
1019 
1020 
1022 
1023  template <int dim>
1024  double
1025  ExpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1026  {
1027  switch (dim)
1028  {
1029  case 1:
1030  return std::exp(p(0));
1031  case 2:
1032  return std::exp(p(0)) * std::exp(p(1));
1033  case 3:
1034  return std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1035  default:
1036  Assert(false, ExcNotImplemented());
1037  }
1038  return 0.;
1039  }
1040 
1041  template <int dim>
1042  void
1043  ExpFunction<dim>::value_list(const std::vector<Point<dim>> &points,
1044  std::vector<double> & values,
1045  const unsigned int) const
1046  {
1047  Assert(values.size() == points.size(),
1048  ExcDimensionMismatch(values.size(), points.size()));
1049 
1050  for (unsigned int i = 0; i < points.size(); ++i)
1051  {
1052  const Point<dim> &p = points[i];
1053  switch (dim)
1054  {
1055  case 1:
1056  values[i] = std::exp(p(0));
1057  break;
1058  case 2:
1059  values[i] = std::exp(p(0)) * std::exp(p(1));
1060  break;
1061  case 3:
1062  values[i] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1063  break;
1064  default:
1065  Assert(false, ExcNotImplemented());
1066  }
1067  }
1068  }
1069 
1070  template <int dim>
1071  double
1072  ExpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1073  {
1074  switch (dim)
1075  {
1076  case 1:
1077  return std::exp(p(0));
1078  case 2:
1079  return 2 * std::exp(p(0)) * std::exp(p(1));
1080  case 3:
1081  return 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1082  default:
1083  Assert(false, ExcNotImplemented());
1084  }
1085  return 0.;
1086  }
1087 
1088  template <int dim>
1089  void
1090  ExpFunction<dim>::laplacian_list(const std::vector<Point<dim>> &points,
1091  std::vector<double> & values,
1092  const unsigned int) const
1093  {
1094  Assert(values.size() == points.size(),
1095  ExcDimensionMismatch(values.size(), points.size()));
1096 
1097  for (unsigned int i = 0; i < points.size(); ++i)
1098  {
1099  const Point<dim> &p = points[i];
1100  switch (dim)
1101  {
1102  case 1:
1103  values[i] = std::exp(p(0));
1104  break;
1105  case 2:
1106  values[i] = 2 * std::exp(p(0)) * std::exp(p(1));
1107  break;
1108  case 3:
1109  values[i] = 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1110  break;
1111  default:
1112  Assert(false, ExcNotImplemented());
1113  }
1114  }
1115  }
1116 
1117  template <int dim>
1118  Tensor<1, dim>
1119  ExpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1120  {
1121  Tensor<1, dim> result;
1122  switch (dim)
1123  {
1124  case 1:
1125  result[0] = std::exp(p(0));
1126  break;
1127  case 2:
1128  result[0] = std::exp(p(0)) * std::exp(p(1));
1129  result[1] = result[0];
1130  break;
1131  case 3:
1132  result[0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1133  result[1] = result[0];
1134  result[2] = result[0];
1135  break;
1136  default:
1137  Assert(false, ExcNotImplemented());
1138  }
1139  return result;
1140  }
1141 
1142  template <int dim>
1143  void
1144  ExpFunction<dim>::gradient_list(const std::vector<Point<dim>> &points,
1145  std::vector<Tensor<1, dim>> & gradients,
1146  const unsigned int) const
1147  {
1148  Assert(gradients.size() == points.size(),
1149  ExcDimensionMismatch(gradients.size(), points.size()));
1150 
1151  for (unsigned int i = 0; i < points.size(); ++i)
1152  {
1153  const Point<dim> &p = points[i];
1154  switch (dim)
1155  {
1156  case 1:
1157  gradients[i][0] = std::exp(p(0));
1158  break;
1159  case 2:
1160  gradients[i][0] = std::exp(p(0)) * std::exp(p(1));
1161  gradients[i][1] = gradients[i][0];
1162  break;
1163  case 3:
1164  gradients[i][0] =
1165  std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
1166  gradients[i][1] = gradients[i][0];
1167  gradients[i][2] = gradients[i][0];
1168  break;
1169  default:
1170  Assert(false, ExcNotImplemented());
1171  }
1172  }
1173  }
1174 
1176 
1177 
1178  double
1179  LSingularityFunction::value(const Point<2> &p, const unsigned int) const
1180  {
1181  double x = p(0);
1182  double y = p(1);
1183 
1184  if ((x >= 0) && (y >= 0))
1185  return 0.;
1186 
1187  double phi = std::atan2(y, -x) + numbers::PI;
1188  double r2 = x * x + y * y;
1189 
1190  return std::pow(r2, 1. / 3.) * std::sin(2. / 3. * phi);
1191  }
1192 
1193 
1194  void
1195  LSingularityFunction::value_list(const std::vector<Point<2>> &points,
1196  std::vector<double> & values,
1197  const unsigned int) const
1198  {
1199  Assert(values.size() == points.size(),
1200  ExcDimensionMismatch(values.size(), points.size()));
1201 
1202  for (unsigned int i = 0; i < points.size(); ++i)
1203  {
1204  double x = points[i](0);
1205  double y = points[i](1);
1206 
1207  if ((x >= 0) && (y >= 0))
1208  values[i] = 0.;
1209  else
1210  {
1211  double phi = std::atan2(y, -x) + numbers::PI;
1212  double r2 = x * x + y * y;
1213 
1214  values[i] = std::pow(r2, 1. / 3.) * std::sin(2. / 3. * phi);
1215  }
1216  }
1217  }
1218 
1219 
1220  void
1221  LSingularityFunction::vector_value_list(
1222  const std::vector<Point<2>> &points,
1223  std::vector<Vector<double>> &values) const
1224  {
1225  Assert(values.size() == points.size(),
1226  ExcDimensionMismatch(values.size(), points.size()));
1227 
1228  for (unsigned int i = 0; i < points.size(); ++i)
1229  {
1230  Assert(values[i].size() == 1,
1231  ExcDimensionMismatch(values[i].size(), 1));
1232  double x = points[i](0);
1233  double y = points[i](1);
1234 
1235  if ((x >= 0) && (y >= 0))
1236  values[i](0) = 0.;
1237  else
1238  {
1239  double phi = std::atan2(y, -x) + numbers::PI;
1240  double r2 = x * x + y * y;
1241 
1242  values[i](0) = std::pow(r2, 1. / 3.) * std::sin(2. / 3. * phi);
1243  }
1244  }
1245  }
1246 
1247 
1248  double
1249  LSingularityFunction::laplacian(const Point<2> &, const unsigned int) const
1250  {
1251  return 0.;
1252  }
1253 
1254 
1255  void
1256  LSingularityFunction::laplacian_list(const std::vector<Point<2>> &points,
1257  std::vector<double> & values,
1258  const unsigned int) const
1259  {
1260  Assert(values.size() == points.size(),
1261  ExcDimensionMismatch(values.size(), points.size()));
1262 
1263  for (unsigned int i = 0; i < points.size(); ++i)
1264  values[i] = 0.;
1265  }
1266 
1267 
1268  Tensor<1, 2>
1269  LSingularityFunction::gradient(const Point<2> &p, const unsigned int) const
1270  {
1271  double x = p(0);
1272  double y = p(1);
1273  double phi = std::atan2(y, -x) + numbers::PI;
1274  double r43 = std::pow(x * x + y * y, 2. / 3.);
1275 
1276  Tensor<1, 2> result;
1277  result[0] = 2. / 3. *
1278  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) /
1279  r43;
1280  result[1] = 2. / 3. *
1281  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) /
1282  r43;
1283  return result;
1284  }
1285 
1286 
1287  void
1288  LSingularityFunction::gradient_list(const std::vector<Point<2>> &points,
1289  std::vector<Tensor<1, 2>> & gradients,
1290  const unsigned int) const
1291  {
1292  Assert(gradients.size() == points.size(),
1293  ExcDimensionMismatch(gradients.size(), points.size()));
1294 
1295  for (unsigned int i = 0; i < points.size(); ++i)
1296  {
1297  const Point<2> &p = points[i];
1298  double x = p(0);
1299  double y = p(1);
1300  double phi = std::atan2(y, -x) + numbers::PI;
1301  double r43 = std::pow(x * x + y * y, 2. / 3.);
1302 
1303  gradients[i][0] =
1304  2. / 3. *
1305  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1306  gradients[i][1] =
1307  2. / 3. *
1308  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1309  }
1310  }
1311 
1312 
1313  void
1314  LSingularityFunction::vector_gradient_list(
1315  const std::vector<Point<2>> & points,
1316  std::vector<std::vector<Tensor<1, 2>>> &gradients) const
1317  {
1318  Assert(gradients.size() == points.size(),
1319  ExcDimensionMismatch(gradients.size(), points.size()));
1320 
1321  for (unsigned int i = 0; i < points.size(); ++i)
1322  {
1323  Assert(gradients[i].size() == 1,
1324  ExcDimensionMismatch(gradients[i].size(), 1));
1325  const Point<2> &p = points[i];
1326  double x = p(0);
1327  double y = p(1);
1328  double phi = std::atan2(y, -x) + numbers::PI;
1329  double r43 = std::pow(x * x + y * y, 2. / 3.);
1330 
1331  gradients[i][0][0] =
1332  2. / 3. *
1333  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1334  gradients[i][0][1] =
1335  2. / 3. *
1336  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1337  }
1338  }
1339 
1341 
1342  LSingularityGradFunction::LSingularityGradFunction()
1343  : Function<2>(2)
1344  {}
1345 
1346 
1347  double
1348  LSingularityGradFunction::value(const Point<2> &p, const unsigned int d) const
1349  {
1350  AssertIndexRange(d, 2);
1351 
1352  const double x = p(0);
1353  const double y = p(1);
1354  const double phi = std::atan2(y, -x) + numbers::PI;
1355  const double r43 = std::pow(x * x + y * y, 2. / 3.);
1356 
1357  return 2. / 3. *
1358  (std::sin(2. / 3. * phi) * p(d) +
1359  (d == 0 ? (std::cos(2. / 3. * phi) * y) :
1360  (-std::cos(2. / 3. * phi) * x))) /
1361  r43;
1362  }
1363 
1364 
1365  void
1366  LSingularityGradFunction::value_list(const std::vector<Point<2>> &points,
1367  std::vector<double> & values,
1368  const unsigned int d) const
1369  {
1370  AssertIndexRange(d, 2);
1371  AssertDimension(values.size(), points.size());
1372 
1373  for (unsigned int i = 0; i < points.size(); ++i)
1374  {
1375  const Point<2> &p = points[i];
1376  const double x = p(0);
1377  const double y = p(1);
1378  const double phi = std::atan2(y, -x) + numbers::PI;
1379  const double r43 = std::pow(x * x + y * y, 2. / 3.);
1380 
1381  values[i] = 2. / 3. *
1382  (std::sin(2. / 3. * phi) * p(d) +
1383  (d == 0 ? (std::cos(2. / 3. * phi) * y) :
1384  (-std::cos(2. / 3. * phi) * x))) /
1385  r43;
1386  }
1387  }
1388 
1389 
1390  void
1391  LSingularityGradFunction::vector_value_list(
1392  const std::vector<Point<2>> &points,
1393  std::vector<Vector<double>> &values) const
1394  {
1395  Assert(values.size() == points.size(),
1396  ExcDimensionMismatch(values.size(), points.size()));
1397 
1398  for (unsigned int i = 0; i < points.size(); ++i)
1399  {
1400  AssertDimension(values[i].size(), 2);
1401  const Point<2> &p = points[i];
1402  const double x = p(0);
1403  const double y = p(1);
1404  const double phi = std::atan2(y, -x) + numbers::PI;
1405  const double r43 = std::pow(x * x + y * y, 2. / 3.);
1406 
1407  values[i](0) =
1408  2. / 3. *
1409  (std::sin(2. / 3. * phi) * x + std::cos(2. / 3. * phi) * y) / r43;
1410  values[i](1) =
1411  2. / 3. *
1412  (std::sin(2. / 3. * phi) * y - std::cos(2. / 3. * phi) * x) / r43;
1413  }
1414  }
1415 
1416 
1417  double
1418  LSingularityGradFunction::laplacian(const Point<2> &,
1419  const unsigned int) const
1420  {
1421  return 0.;
1422  }
1423 
1424 
1425  void
1426  LSingularityGradFunction::laplacian_list(const std::vector<Point<2>> &points,
1427  std::vector<double> & values,
1428  const unsigned int) const
1429  {
1430  Assert(values.size() == points.size(),
1431  ExcDimensionMismatch(values.size(), points.size()));
1432 
1433  for (unsigned int i = 0; i < points.size(); ++i)
1434  values[i] = 0.;
1435  }
1436 
1437 
1438 
1439  Tensor<1, 2>
1440  LSingularityGradFunction::gradient(const Point<2> & /*p*/,
1441  const unsigned int /*component*/) const
1442  {
1443  Assert(false, ExcNotImplemented());
1444  return Tensor<1, 2>();
1445  }
1446 
1447 
1448  void
1449  LSingularityGradFunction::gradient_list(
1450  const std::vector<Point<2>> & /*points*/,
1451  std::vector<Tensor<1, 2>> & /*gradients*/,
1452  const unsigned int /*component*/) const
1453  {
1454  Assert(false, ExcNotImplemented());
1455  }
1456 
1457 
1458  void
1459  LSingularityGradFunction::vector_gradient_list(
1460  const std::vector<Point<2>> & /*points*/,
1461  std::vector<std::vector<Tensor<1, 2>>> & /*gradients*/) const
1462  {
1463  Assert(false, ExcNotImplemented());
1464  }
1465 
1467 
1468  template <int dim>
1469  double
1470  SlitSingularityFunction<dim>::value(const Point<dim> &p,
1471  const unsigned int) const
1472  {
1473  double x = p(0);
1474  double y = p(1);
1475 
1476  double phi = std::atan2(x, y) + numbers::PI;
1477  double r2 = x * x + y * y;
1478 
1479  return std::pow(r2, .25) * std::sin(.5 * phi);
1480  }
1481 
1482 
1483  template <int dim>
1484  void
1485  SlitSingularityFunction<dim>::value_list(
1486  const std::vector<Point<dim>> &points,
1487  std::vector<double> & values,
1488  const unsigned int) const
1489  {
1490  Assert(values.size() == points.size(),
1491  ExcDimensionMismatch(values.size(), points.size()));
1492 
1493  for (unsigned int i = 0; i < points.size(); ++i)
1494  {
1495  double x = points[i](0);
1496  double y = points[i](1);
1497 
1498  double phi = std::atan2(x, y) + numbers::PI;
1499  double r2 = x * x + y * y;
1500 
1501  values[i] = std::pow(r2, .25) * std::sin(.5 * phi);
1502  }
1503  }
1504 
1505 
1506  template <int dim>
1507  void
1508  SlitSingularityFunction<dim>::vector_value_list(
1509  const std::vector<Point<dim>> &points,
1510  std::vector<Vector<double>> & values) const
1511  {
1512  Assert(values.size() == points.size(),
1513  ExcDimensionMismatch(values.size(), points.size()));
1514 
1515  for (unsigned int i = 0; i < points.size(); ++i)
1516  {
1517  Assert(values[i].size() == 1,
1518  ExcDimensionMismatch(values[i].size(), 1));
1519 
1520  double x = points[i](0);
1521  double y = points[i](1);
1522 
1523  double phi = std::atan2(x, y) + numbers::PI;
1524  double r2 = x * x + y * y;
1525 
1526  values[i](0) = std::pow(r2, .25) * std::sin(.5 * phi);
1527  }
1528  }
1529 
1530 
1531  template <int dim>
1532  double
1533  SlitSingularityFunction<dim>::laplacian(const Point<dim> &,
1534  const unsigned int) const
1535  {
1536  return 0.;
1537  }
1538 
1539 
1540  template <int dim>
1541  void
1542  SlitSingularityFunction<dim>::laplacian_list(
1543  const std::vector<Point<dim>> &points,
1544  std::vector<double> & values,
1545  const unsigned int) const
1546  {
1547  Assert(values.size() == points.size(),
1548  ExcDimensionMismatch(values.size(), points.size()));
1549 
1550  for (unsigned int i = 0; i < points.size(); ++i)
1551  values[i] = 0.;
1552  }
1553 
1554 
1555  template <int dim>
1556  Tensor<1, dim>
1557  SlitSingularityFunction<dim>::gradient(const Point<dim> &p,
1558  const unsigned int) const
1559  {
1560  double x = p(0);
1561  double y = p(1);
1562  double phi = std::atan2(x, y) + numbers::PI;
1563  double r64 = std::pow(x * x + y * y, 3. / 4.);
1564 
1565  Tensor<1, dim> result;
1566  result[0] = 1. / 2. *
1567  (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) /
1568  r64;
1569  result[1] = 1. / 2. *
1570  (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) /
1571  r64;
1572  return result;
1573  }
1574 
1575 
1576  template <int dim>
1577  void
1578  SlitSingularityFunction<dim>::gradient_list(
1579  const std::vector<Point<dim>> &points,
1580  std::vector<Tensor<1, dim>> & gradients,
1581  const unsigned int) const
1582  {
1583  Assert(gradients.size() == points.size(),
1584  ExcDimensionMismatch(gradients.size(), points.size()));
1585 
1586  for (unsigned int i = 0; i < points.size(); ++i)
1587  {
1588  const Point<dim> &p = points[i];
1589  double x = p(0);
1590  double y = p(1);
1591  double phi = std::atan2(x, y) + numbers::PI;
1592  double r64 = std::pow(x * x + y * y, 3. / 4.);
1593 
1594  gradients[i][0] =
1595  1. / 2. *
1596  (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) / r64;
1597  gradients[i][1] =
1598  1. / 2. *
1599  (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) / r64;
1600  for (unsigned int d = 2; d < dim; ++d)
1601  gradients[i][d] = 0.;
1602  }
1603  }
1604 
1605  template <int dim>
1606  void
1607  SlitSingularityFunction<dim>::vector_gradient_list(
1608  const std::vector<Point<dim>> & points,
1609  std::vector<std::vector<Tensor<1, dim>>> &gradients) const
1610  {
1611  Assert(gradients.size() == points.size(),
1612  ExcDimensionMismatch(gradients.size(), points.size()));
1613 
1614  for (unsigned int i = 0; i < points.size(); ++i)
1615  {
1616  Assert(gradients[i].size() == 1,
1617  ExcDimensionMismatch(gradients[i].size(), 1));
1618 
1619  const Point<dim> &p = points[i];
1620  double x = p(0);
1621  double y = p(1);
1622  double phi = std::atan2(x, y) + numbers::PI;
1623  double r64 = std::pow(x * x + y * y, 3. / 4.);
1624 
1625  gradients[i][0][0] =
1626  1. / 2. *
1627  (std::sin(1. / 2. * phi) * x + std::cos(1. / 2. * phi) * y) / r64;
1628  gradients[i][0][1] =
1629  1. / 2. *
1630  (std::sin(1. / 2. * phi) * y - std::cos(1. / 2. * phi) * x) / r64;
1631  for (unsigned int d = 2; d < dim; ++d)
1632  gradients[i][0][d] = 0.;
1633  }
1634  }
1635 
1637 
1638 
1639  double
1640  SlitHyperSingularityFunction::value(const Point<2> &p,
1641  const unsigned int) const
1642  {
1643  double x = p(0);
1644  double y = p(1);
1645 
1646  double phi = std::atan2(x, y) + numbers::PI;
1647  double r2 = x * x + y * y;
1648 
1649  return std::pow(r2, .125) * std::sin(.25 * phi);
1650  }
1651 
1652 
1653  void
1654  SlitHyperSingularityFunction::value_list(const std::vector<Point<2>> &points,
1655  std::vector<double> & values,
1656  const unsigned int) const
1657  {
1658  Assert(values.size() == points.size(),
1659  ExcDimensionMismatch(values.size(), points.size()));
1660 
1661  for (unsigned int i = 0; i < points.size(); ++i)
1662  {
1663  double x = points[i](0);
1664  double y = points[i](1);
1665 
1666  double phi = std::atan2(x, y) + numbers::PI;
1667  double r2 = x * x + y * y;
1668 
1669  values[i] = std::pow(r2, .125) * std::sin(.25 * phi);
1670  }
1671  }
1672 
1673 
1674  void
1675  SlitHyperSingularityFunction::vector_value_list(
1676  const std::vector<Point<2>> &points,
1677  std::vector<Vector<double>> &values) const
1678  {
1679  Assert(values.size() == points.size(),
1680  ExcDimensionMismatch(values.size(), points.size()));
1681 
1682  for (unsigned int i = 0; i < points.size(); ++i)
1683  {
1684  Assert(values[i].size() == 1,
1685  ExcDimensionMismatch(values[i].size(), 1));
1686 
1687  double x = points[i](0);
1688  double y = points[i](1);
1689 
1690  double phi = std::atan2(x, y) + numbers::PI;
1691  double r2 = x * x + y * y;
1692 
1693  values[i](0) = std::pow(r2, .125) * std::sin(.25 * phi);
1694  }
1695  }
1696 
1697 
1698  double
1699  SlitHyperSingularityFunction::laplacian(const Point<2> &,
1700  const unsigned int) const
1701  {
1702  return 0.;
1703  }
1704 
1705 
1706  void
1707  SlitHyperSingularityFunction::laplacian_list(
1708  const std::vector<Point<2>> &points,
1709  std::vector<double> & values,
1710  const unsigned int) const
1711  {
1712  Assert(values.size() == points.size(),
1713  ExcDimensionMismatch(values.size(), points.size()));
1714 
1715  for (unsigned int i = 0; i < points.size(); ++i)
1716  values[i] = 0.;
1717  }
1718 
1719 
1720  Tensor<1, 2>
1721  SlitHyperSingularityFunction::gradient(const Point<2> &p,
1722  const unsigned int) const
1723  {
1724  double x = p(0);
1725  double y = p(1);
1726  double phi = std::atan2(x, y) + numbers::PI;
1727  double r78 = std::pow(x * x + y * y, 7. / 8.);
1728 
1729 
1730  Tensor<1, 2> result;
1731  result[0] = 1. / 4. *
1732  (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) /
1733  r78;
1734  result[1] = 1. / 4. *
1735  (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) /
1736  r78;
1737  return result;
1738  }
1739 
1740 
1741  void
1742  SlitHyperSingularityFunction::gradient_list(
1743  const std::vector<Point<2>> &points,
1744  std::vector<Tensor<1, 2>> & gradients,
1745  const unsigned int) const
1746  {
1747  Assert(gradients.size() == points.size(),
1748  ExcDimensionMismatch(gradients.size(), points.size()));
1749 
1750  for (unsigned int i = 0; i < points.size(); ++i)
1751  {
1752  const Point<2> &p = points[i];
1753  double x = p(0);
1754  double y = p(1);
1755  double phi = std::atan2(x, y) + numbers::PI;
1756  double r78 = std::pow(x * x + y * y, 7. / 8.);
1757 
1758  gradients[i][0] =
1759  1. / 4. *
1760  (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) / r78;
1761  gradients[i][1] =
1762  1. / 4. *
1763  (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) / r78;
1764  }
1765  }
1766 
1767 
1768  void
1769  SlitHyperSingularityFunction::vector_gradient_list(
1770  const std::vector<Point<2>> & points,
1771  std::vector<std::vector<Tensor<1, 2>>> &gradients) const
1772  {
1773  Assert(gradients.size() == points.size(),
1774  ExcDimensionMismatch(gradients.size(), points.size()));
1775 
1776  for (unsigned int i = 0; i < points.size(); ++i)
1777  {
1778  Assert(gradients[i].size() == 1,
1779  ExcDimensionMismatch(gradients[i].size(), 1));
1780 
1781  const Point<2> &p = points[i];
1782  double x = p(0);
1783  double y = p(1);
1784  double phi = std::atan2(x, y) + numbers::PI;
1785  double r78 = std::pow(x * x + y * y, 7. / 8.);
1786 
1787  gradients[i][0][0] =
1788  1. / 4. *
1789  (std::sin(1. / 4. * phi) * x + std::cos(1. / 4. * phi) * y) / r78;
1790  gradients[i][0][1] =
1791  1. / 4. *
1792  (std::sin(1. / 4. * phi) * y - std::cos(1. / 4. * phi) * x) / r78;
1793  }
1794  }
1795 
1797 
1798  template <int dim>
1799  JumpFunction<dim>::JumpFunction(const Point<dim> &direction,
1800  const double steepness)
1801  : direction(direction)
1802  , steepness(steepness)
1803  {
1804  switch (dim)
1805  {
1806  case 1:
1807  angle = 0;
1808  break;
1809  case 2:
1810  angle = std::atan2(direction(0), direction(1));
1811  break;
1812  case 3:
1813  Assert(false, ExcNotImplemented());
1814  }
1815  sine = std::sin(angle);
1816  cosine = std::cos(angle);
1817  }
1818 
1819 
1820 
1821  template <int dim>
1822  double
1823  JumpFunction<dim>::value(const Point<dim> &p, const unsigned int) const
1824  {
1825  double x = steepness * (-cosine * p(0) + sine * p(1));
1826  return -std::atan(x);
1827  }
1828 
1829 
1830 
1831  template <int dim>
1832  void
1833  JumpFunction<dim>::value_list(const std::vector<Point<dim>> &p,
1834  std::vector<double> & values,
1835  const unsigned int) const
1836  {
1837  Assert(values.size() == p.size(),
1838  ExcDimensionMismatch(values.size(), p.size()));
1839 
1840  for (unsigned int i = 0; i < p.size(); ++i)
1841  {
1842  double x = steepness * (-cosine * p[i](0) + sine * p[i](1));
1843  values[i] = -std::atan(x);
1844  }
1845  }
1846 
1847 
1848  template <int dim>
1849  double
1850  JumpFunction<dim>::laplacian(const Point<dim> &p, const unsigned int) const
1851  {
1852  double x = steepness * (-cosine * p(0) + sine * p(1));
1853  double r = 1 + x * x;
1854  return 2 * steepness * steepness * x / (r * r);
1855  }
1856 
1857 
1858  template <int dim>
1859  void
1860  JumpFunction<dim>::laplacian_list(const std::vector<Point<dim>> &p,
1861  std::vector<double> & values,
1862  const unsigned int) const
1863  {
1864  Assert(values.size() == p.size(),
1865  ExcDimensionMismatch(values.size(), p.size()));
1866 
1867  double f = 2 * steepness * steepness;
1868 
1869  for (unsigned int i = 0; i < p.size(); ++i)
1870  {
1871  double x = steepness * (-cosine * p[i](0) + sine * p[i](1));
1872  double r = 1 + x * x;
1873  values[i] = f * x / (r * r);
1874  }
1875  }
1876 
1877 
1878 
1879  template <int dim>
1880  Tensor<1, dim>
1881  JumpFunction<dim>::gradient(const Point<dim> &p, const unsigned int) const
1882  {
1883  double x = steepness * (-cosine * p(0) + sine * p(1));
1884  double r = -steepness * (1 + x * x);
1885  Tensor<1, dim> erg;
1886  erg[0] = cosine * r;
1887  erg[1] = sine * r;
1888  return erg;
1889  }
1890 
1891 
1892 
1893  template <int dim>
1894  void
1895  JumpFunction<dim>::gradient_list(const std::vector<Point<dim>> &p,
1896  std::vector<Tensor<1, dim>> & gradients,
1897  const unsigned int) const
1898  {
1899  Assert(gradients.size() == p.size(),
1900  ExcDimensionMismatch(gradients.size(), p.size()));
1901 
1902  for (unsigned int i = 0; i < p.size(); ++i)
1903  {
1904  double x = steepness * (cosine * p[i](0) + sine * p[i](1));
1905  double r = -steepness * (1 + x * x);
1906  gradients[i][0] = cosine * r;
1907  gradients[i][1] = sine * r;
1908  }
1909  }
1910 
1911 
1912 
1913  template <int dim>
1914  std::size_t
1915  JumpFunction<dim>::memory_consumption() const
1916  {
1917  // only simple data elements, so
1918  // use sizeof operator
1919  return sizeof(*this);
1920  }
1921 
1922 
1923 
1924  /* ---------------------- FourierCosineFunction ----------------------- */
1925 
1926 
1927  template <int dim>
1928  FourierCosineFunction<dim>::FourierCosineFunction(
1929  const Tensor<1, dim> &fourier_coefficients)
1930  : Function<dim>(1)
1931  , fourier_coefficients(fourier_coefficients)
1932  {}
1933 
1934 
1935 
1936  template <int dim>
1937  double
1938  FourierCosineFunction<dim>::value(const Point<dim> & p,
1939  const unsigned int component) const
1940  {
1941  (void)component;
1942  Assert(component == 0, ExcIndexRange(component, 0, 1));
1943  return std::cos(fourier_coefficients * p);
1944  }
1945 
1946 
1947 
1948  template <int dim>
1949  Tensor<1, dim>
1950  FourierCosineFunction<dim>::gradient(const Point<dim> & p,
1951  const unsigned int component) const
1952  {
1953  (void)component;
1954  Assert(component == 0, ExcIndexRange(component, 0, 1));
1955  return -fourier_coefficients * std::sin(fourier_coefficients * p);
1956  }
1957 
1958 
1959 
1960  template <int dim>
1961  double
1962  FourierCosineFunction<dim>::laplacian(const Point<dim> & p,
1963  const unsigned int component) const
1964  {
1965  (void)component;
1966  Assert(component == 0, ExcIndexRange(component, 0, 1));
1967  return (fourier_coefficients * fourier_coefficients) *
1968  (-std::cos(fourier_coefficients * p));
1969  }
1970 
1971 
1972 
1973  /* ---------------------- FourierSineFunction ----------------------- */
1974 
1975 
1976 
1977  template <int dim>
1978  FourierSineFunction<dim>::FourierSineFunction(
1979  const Tensor<1, dim> &fourier_coefficients)
1980  : Function<dim>(1)
1981  , fourier_coefficients(fourier_coefficients)
1982  {}
1983 
1984 
1985 
1986  template <int dim>
1987  double
1988  FourierSineFunction<dim>::value(const Point<dim> & p,
1989  const unsigned int component) const
1990  {
1991  (void)component;
1992  Assert(component == 0, ExcIndexRange(component, 0, 1));
1993  return std::sin(fourier_coefficients * p);
1994  }
1995 
1996 
1997 
1998  template <int dim>
1999  Tensor<1, dim>
2000  FourierSineFunction<dim>::gradient(const Point<dim> & p,
2001  const unsigned int component) const
2002  {
2003  (void)component;
2004  Assert(component == 0, ExcIndexRange(component, 0, 1));
2005  return fourier_coefficients * std::cos(fourier_coefficients * p);
2006  }
2007 
2008 
2009 
2010  template <int dim>
2011  double
2012  FourierSineFunction<dim>::laplacian(const Point<dim> & p,
2013  const unsigned int component) const
2014  {
2015  (void)component;
2016  Assert(component == 0, ExcIndexRange(component, 0, 1));
2017  return (fourier_coefficients * fourier_coefficients) *
2018  (-std::sin(fourier_coefficients * p));
2019  }
2020 
2021 
2022 
2023  /* ---------------------- FourierSineSum ----------------------- */
2024 
2025 
2026 
2027  template <int dim>
2028  FourierSineSum<dim>::FourierSineSum(
2029  const std::vector<Point<dim>> &fourier_coefficients,
2030  const std::vector<double> & weights)
2031  : Function<dim>(1)
2032  , fourier_coefficients(fourier_coefficients)
2033  , weights(weights)
2034  {
2035  Assert(fourier_coefficients.size() > 0, ExcZero());
2036  Assert(fourier_coefficients.size() == weights.size(),
2037  ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2038  }
2039 
2040 
2041 
2042  template <int dim>
2043  double
2044  FourierSineSum<dim>::value(const Point<dim> & p,
2045  const unsigned int component) const
2046  {
2047  (void)component;
2048  Assert(component == 0, ExcIndexRange(component, 0, 1));
2049 
2050  const unsigned int n = weights.size();
2051  double sum = 0;
2052  for (unsigned int s = 0; s < n; ++s)
2053  sum += weights[s] * std::sin(fourier_coefficients[s] * p);
2054 
2055  return sum;
2056  }
2057 
2058 
2059 
2060  template <int dim>
2061  Tensor<1, dim>
2062  FourierSineSum<dim>::gradient(const Point<dim> & p,
2063  const unsigned int component) const
2064  {
2065  (void)component;
2066  Assert(component == 0, ExcIndexRange(component, 0, 1));
2067 
2068  const unsigned int n = weights.size();
2069  Tensor<1, dim> sum;
2070  for (unsigned int s = 0; s < n; ++s)
2071  sum += fourier_coefficients[s] * std::cos(fourier_coefficients[s] * p);
2072 
2073  return sum;
2074  }
2075 
2076 
2077 
2078  template <int dim>
2079  double
2080  FourierSineSum<dim>::laplacian(const Point<dim> & p,
2081  const unsigned int component) const
2082  {
2083  (void)component;
2084  Assert(component == 0, ExcIndexRange(component, 0, 1));
2085 
2086  const unsigned int n = weights.size();
2087  double sum = 0;
2088  for (unsigned int s = 0; s < n; ++s)
2089  sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2090  std::sin(fourier_coefficients[s] * p);
2091 
2092  return sum;
2093  }
2094 
2095 
2096 
2097  /* ---------------------- FourierCosineSum ----------------------- */
2098 
2099 
2100 
2101  template <int dim>
2102  FourierCosineSum<dim>::FourierCosineSum(
2103  const std::vector<Point<dim>> &fourier_coefficients,
2104  const std::vector<double> & weights)
2105  : Function<dim>(1)
2106  , fourier_coefficients(fourier_coefficients)
2107  , weights(weights)
2108  {
2109  Assert(fourier_coefficients.size() > 0, ExcZero());
2110  Assert(fourier_coefficients.size() == weights.size(),
2111  ExcDimensionMismatch(fourier_coefficients.size(), weights.size()));
2112  }
2113 
2114 
2115 
2116  template <int dim>
2117  double
2118  FourierCosineSum<dim>::value(const Point<dim> & p,
2119  const unsigned int component) const
2120  {
2121  (void)component;
2122  Assert(component == 0, ExcIndexRange(component, 0, 1));
2123 
2124  const unsigned int n = weights.size();
2125  double sum = 0;
2126  for (unsigned int s = 0; s < n; ++s)
2127  sum += weights[s] * std::cos(fourier_coefficients[s] * p);
2128 
2129  return sum;
2130  }
2131 
2132 
2133 
2134  template <int dim>
2135  Tensor<1, dim>
2136  FourierCosineSum<dim>::gradient(const Point<dim> & p,
2137  const unsigned int component) const
2138  {
2139  (void)component;
2140  Assert(component == 0, ExcIndexRange(component, 0, 1));
2141 
2142  const unsigned int n = weights.size();
2143  Tensor<1, dim> sum;
2144  for (unsigned int s = 0; s < n; ++s)
2145  sum -= fourier_coefficients[s] * std::sin(fourier_coefficients[s] * p);
2146 
2147  return sum;
2148  }
2149 
2150 
2151 
2152  template <int dim>
2153  double
2154  FourierCosineSum<dim>::laplacian(const Point<dim> & p,
2155  const unsigned int component) const
2156  {
2157  (void)component;
2158  Assert(component == 0, ExcIndexRange(component, 0, 1));
2159 
2160  const unsigned int n = weights.size();
2161  double sum = 0;
2162  for (unsigned int s = 0; s < n; ++s)
2163  sum -= (fourier_coefficients[s] * fourier_coefficients[s]) *
2164  std::cos(fourier_coefficients[s] * p);
2165 
2166  return sum;
2167  }
2168 
2169 
2170 
2171  /* ---------------------- Monomial ----------------------- */
2172 
2173 
2174 
2175  template <int dim>
2176  Monomial<dim>::Monomial(const Tensor<1, dim> &exponents,
2177  const unsigned int n_components)
2178  : Function<dim>(n_components)
2179  , exponents(exponents)
2180  {}
2181 
2182 
2183 
2184  template <int dim>
2185  double
2186  Monomial<dim>::value(const Point<dim> &p, const unsigned int component) const
2187  {
2188  (void)component;
2189  Assert(component < this->n_components,
2190  ExcIndexRange(component, 0, this->n_components));
2191 
2192  double prod = 1;
2193  for (unsigned int s = 0; s < dim; ++s)
2194  {
2195  if (p[s] < 0)
2196  Assert(std::floor(exponents[s]) == exponents[s],
2197  ExcMessage("Exponentiation of a negative base number with "
2198  "a real exponent can't be performed."));
2199  prod *= std::pow(p[s], exponents[s]);
2200  }
2201  return prod;
2202  }
2203 
2204 
2205 
2206  template <int dim>
2207  void
2208  Monomial<dim>::vector_value(const Point<dim> &p, Vector<double> &values) const
2209  {
2210  Assert(values.size() == this->n_components,
2211  ExcDimensionMismatch(values.size(), this->n_components));
2212 
2213  for (unsigned int i = 0; i < values.size(); ++i)
2214  values(i) = Monomial<dim>::value(p, i);
2215  }
2216 
2217 
2218 
2219  template <int dim>
2220  Tensor<1, dim>
2221  Monomial<dim>::gradient(const Point<dim> & p,
2222  const unsigned int component) const
2223  {
2224  (void)component;
2225  Assert(component == 0, ExcIndexRange(component, 0, 1));
2226 
2227  Tensor<1, dim> r;
2228  for (unsigned int d = 0; d < dim; ++d)
2229  {
2230  double prod = 1;
2231  for (unsigned int s = 0; s < dim; ++s)
2232  {
2233  if ((s == d) && (exponents[s] == 0) && (p[s] == 0))
2234  {
2235  prod = 0;
2236  break;
2237  }
2238  else
2239  {
2240  if (p[s] < 0)
2241  Assert(std::floor(exponents[s]) == exponents[s],
2242  ExcMessage(
2243  "Exponentiation of a negative base number with "
2244  "a real exponent can't be performed."));
2245  prod *=
2246  (s == d ? exponents[s] * std::pow(p[s], exponents[s] - 1) :
2247  std::pow(p[s], exponents[s]));
2248  }
2249  }
2250  r[d] = prod;
2251  }
2252 
2253  return r;
2254  }
2255 
2256 
2257 
2258  template <int dim>
2259  void
2260  Monomial<dim>::value_list(const std::vector<Point<dim>> &points,
2261  std::vector<double> & values,
2262  const unsigned int component) const
2263  {
2264  Assert(values.size() == points.size(),
2265  ExcDimensionMismatch(values.size(), points.size()));
2266 
2267  for (unsigned int i = 0; i < points.size(); ++i)
2268  values[i] = Monomial<dim>::value(points[i], component);
2269  }
2270 
2271 
2272  template <int dim>
2273  Bessel1<dim>::Bessel1(const unsigned int order,
2274  const double wave_number,
2275  const Point<dim> center)
2276  : order(order)
2277  , wave_number(wave_number)
2278  , center(center)
2279  {
2280  Assert(wave_number >= 0., ExcMessage("wave_number must be nonnegative!"));
2281  }
2282 
2283  template <int dim>
2284  double
2285  Bessel1<dim>::value(const Point<dim> &p, const unsigned int) const
2286  {
2287  Assert(dim == 2, ExcNotImplemented());
2288  const double r = p.distance(center);
2289  return std_cxx17::cyl_bessel_j(order, r * wave_number);
2290  }
2291 
2292 
2293  template <int dim>
2294  void
2295  Bessel1<dim>::value_list(const std::vector<Point<dim>> &points,
2296  std::vector<double> & values,
2297  const unsigned int) const
2298  {
2299  Assert(dim == 2, ExcNotImplemented());
2300  AssertDimension(points.size(), values.size());
2301  for (unsigned int k = 0; k < points.size(); ++k)
2302  {
2303  const double r = points[k].distance(center);
2304  values[k] = std_cxx17::cyl_bessel_j(order, r * wave_number);
2305  }
2306  }
2307 
2308 
2309  template <int dim>
2310  Tensor<1, dim>
2311  Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
2312  {
2313  Assert(dim == 2, ExcNotImplemented());
2314  const double r = p.distance(center);
2315  const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
2316  const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
2317 
2318  const double dJn =
2319  (order == 0) ?
2320  (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2321  (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2322  std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2323  Tensor<1, dim> result;
2324  result[0] = wave_number * co * dJn;
2325  result[1] = wave_number * si * dJn;
2326  return result;
2327  }
2328 
2329 
2330 
2331  template <int dim>
2332  void
2333  Bessel1<dim>::gradient_list(const std::vector<Point<dim>> &points,
2334  std::vector<Tensor<1, dim>> & gradients,
2335  const unsigned int) const
2336  {
2337  Assert(dim == 2, ExcNotImplemented());
2338  AssertDimension(points.size(), gradients.size());
2339  for (unsigned int k = 0; k < points.size(); ++k)
2340  {
2341  const Point<dim> &p = points[k];
2342  const double r = p.distance(center);
2343  const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
2344  const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
2345 
2346  const double dJn =
2347  (order == 0) ?
2348  (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
2349  (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
2350  std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
2351  Tensor<1, dim> &result = gradients[k];
2352  result[0] = wave_number * co * dJn;
2353  result[1] = wave_number * si * dJn;
2354  }
2355  }
2356 
2357 
2358 
2359  namespace
2360  {
2361  // interpolate a data value from a table where ix denotes
2362  // the (lower) left endpoint of the interval to interpolate
2363  // in, and p_unit denotes the point in unit coordinates to do so.
2364  double
2365  interpolate(const Table<1, double> &data_values,
2366  const TableIndices<1> & ix,
2367  const Point<1> & xi)
2368  {
2369  return ((1 - xi[0]) * data_values[ix[0]] +
2370  xi[0] * data_values[ix[0] + 1]);
2371  }
2372 
2373  double
2374  interpolate(const Table<2, double> &data_values,
2375  const TableIndices<2> & ix,
2376  const Point<2> & p_unit)
2377  {
2378  return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
2379  p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
2380  (1 - p_unit[1]) +
2381  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
2382  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
2383  p_unit[1]);
2384  }
2385 
2386  double
2387  interpolate(const Table<3, double> &data_values,
2388  const TableIndices<3> & ix,
2389  const Point<3> & p_unit)
2390  {
2391  return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
2392  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
2393  (1 - p_unit[1]) +
2394  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
2395  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
2396  p_unit[1]) *
2397  (1 - p_unit[2]) +
2398  (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
2399  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
2400  (1 - p_unit[1]) +
2401  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
2402  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
2403  p_unit[1]) *
2404  p_unit[2]);
2405  }
2406 
2407 
2408  // Interpolate the gradient of a data value from a table where ix
2409  // denotes the lower left endpoint of the interval to interpolate
2410  // in, p_unit denotes the point in unit coordinates, and dx
2411  // denotes the width of the interval in each dimension.
2412  Tensor<1, 1>
2413  gradient_interpolate(const Table<1, double> &data_values,
2414  const TableIndices<1> & ix,
2415  const Point<1> & p_unit,
2416  const Point<1> & dx)
2417  {
2418  (void)p_unit;
2419  Tensor<1, 1> grad;
2420  grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
2421  return grad;
2422  }
2423 
2424 
2425  Tensor<1, 2>
2426  gradient_interpolate(const Table<2, double> &data_values,
2427  const TableIndices<2> & ix,
2428  const Point<2> & p_unit,
2429  const Point<2> & dx)
2430  {
2431  Tensor<1, 2> grad;
2432  double u00 = data_values[ix[0]][ix[1]],
2433  u01 = data_values[ix[0] + 1][ix[1]],
2434  u10 = data_values[ix[0]][ix[1] + 1],
2435  u11 = data_values[ix[0] + 1][ix[1] + 1];
2436 
2437  grad[0] =
2438  ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
2439  grad[1] =
2440  ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
2441  return grad;
2442  }
2443 
2444 
2445  Tensor<1, 3>
2446  gradient_interpolate(const Table<3, double> &data_values,
2447  const TableIndices<3> & ix,
2448  const Point<3> & p_unit,
2449  const Point<3> & dx)
2450  {
2451  Tensor<1, 3> grad;
2452  double u000 = data_values[ix[0]][ix[1]][ix[2]],
2453  u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
2454  u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
2455  u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
2456  u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
2457  u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
2458  u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
2459  u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
2460 
2461  grad[0] =
2462  ((1 - p_unit[2]) *
2463  ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
2464  p_unit[2] *
2465  ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
2466  dx[0];
2467  grad[1] =
2468  ((1 - p_unit[2]) *
2469  ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
2470  p_unit[2] *
2471  ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
2472  dx[1];
2473  grad[2] =
2474  ((1 - p_unit[1]) *
2475  ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
2476  p_unit[1] *
2477  ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
2478  dx[2];
2479 
2480  return grad;
2481  }
2482  } // namespace
2483 
2484  template <int dim>
2485  InterpolatedTensorProductGridData<dim>::InterpolatedTensorProductGridData(
2486  const std::array<std::vector<double>, dim> &coordinate_values,
2487  const Table<dim, double> & data_values)
2488  : coordinate_values(coordinate_values)
2489  , data_values(data_values)
2490  {
2491  for (unsigned int d = 0; d < dim; ++d)
2492  {
2493  Assert(
2494  coordinate_values[d].size() >= 2,
2495  ExcMessage(
2496  "Coordinate arrays must have at least two coordinate values!"));
2497  for (unsigned int i = 0; i < coordinate_values[d].size() - 1; ++i)
2498  Assert(
2499  coordinate_values[d][i] < coordinate_values[d][i + 1],
2500  ExcMessage(
2501  "Coordinate arrays must be sorted in strictly ascending order."));
2502 
2503  Assert(data_values.size()[d] == coordinate_values[d].size(),
2504  ExcMessage(
2505  "Data and coordinate tables do not have the same size."));
2506  }
2507  }
2508 
2509 
2510 
2511  template <int dim>
2512  TableIndices<dim>
2513  InterpolatedTensorProductGridData<dim>::table_index_of_point(
2514  const Point<dim> &p) const
2515  {
2516  // find out where this data point lies, relative to the given
2517  // points. if we run all the way to the end of the range,
2518  // set the indices so that we will simply query the last of the
2519  // intervals, starting at x.size()-2 and going to x.size()-1.
2520  TableIndices<dim> ix;
2521  for (unsigned int d = 0; d < dim; ++d)
2522  {
2523  // get the index of the first element of the coordinate arrays that is
2524  // larger than p[d]
2525  ix[d] = (std::lower_bound(coordinate_values[d].begin(),
2526  coordinate_values[d].end(),
2527  p[d]) -
2528  coordinate_values[d].begin());
2529 
2530  // the one we want is the index of the coordinate to the left, however,
2531  // so decrease it by one (unless we have a point to the left of all, in
2532  // which case we stay where we are; the formulas below are made in a way
2533  // that allow us to extend the function by a constant value)
2534  //
2535  // to make this work, if we got coordinate_values[d].end(), we actually
2536  // have to consider the last box which has index size()-2
2537  if (ix[d] == coordinate_values[d].size())
2538  ix[d] = coordinate_values[d].size() - 2;
2539  else if (ix[d] > 0)
2540  --ix[d];
2541  }
2542 
2543  return ix;
2544  }
2545 
2546 
2547 
2548  template <int dim>
2549  double
2550  InterpolatedTensorProductGridData<dim>::value(
2551  const Point<dim> & p,
2552  const unsigned int component) const
2553  {
2554  (void)component;
2555  Assert(
2556  component == 0,
2557  ExcMessage(
2558  "This is a scalar function object, the component can only be zero."));
2559 
2560  // find the index in the data table of the cell containing the input point
2561  const TableIndices<dim> ix = table_index_of_point(p);
2562 
2563  // now compute the relative point within the interval/rectangle/box
2564  // defined by the point coordinates found above. truncate below and
2565  // above to accommodate points that may lie outside the range
2566  Point<dim> p_unit;
2567  for (unsigned int d = 0; d < dim; ++d)
2568  p_unit[d] = std::max(std::min((p[d] - coordinate_values[d][ix[d]]) /
2569  (coordinate_values[d][ix[d] + 1] -
2570  coordinate_values[d][ix[d]]),
2571  1.),
2572  0.);
2573 
2574  return interpolate(data_values, ix, p_unit);
2575  }
2576 
2577 
2578 
2579  template <int dim>
2580  Tensor<1, dim>
2581  InterpolatedTensorProductGridData<dim>::gradient(
2582  const Point<dim> & p,
2583  const unsigned int component) const
2584  {
2585  (void)component;
2586  Assert(
2587  component == 0,
2588  ExcMessage(
2589  "This is a scalar function object, the component can only be zero."));
2590 
2591  // find out where this data point lies
2592  const TableIndices<dim> ix = table_index_of_point(p);
2593 
2594  Point<dim> dx;
2595  for (unsigned int d = 0; d < dim; ++d)
2596  dx[d] = coordinate_values[d][ix[d] + 1] - coordinate_values[d][ix[d]];
2597 
2598  Point<dim> p_unit;
2599  for (unsigned int d = 0; d < dim; ++d)
2600  p_unit[d] =
2601  std::max(std::min((p[d] - coordinate_values[d][ix[d]]) / dx[d], 1.),
2602  0.0);
2603 
2604  return gradient_interpolate(data_values, ix, p_unit, dx);
2605  }
2606 
2607 
2608 
2609  template <int dim>
2610  InterpolatedUniformGridData<dim>::InterpolatedUniformGridData(
2611  const std::array<std::pair<double, double>, dim> &interval_endpoints,
2612  const std::array<unsigned int, dim> & n_subintervals,
2613  const Table<dim, double> & data_values)
2614  : interval_endpoints(interval_endpoints)
2615  , n_subintervals(n_subintervals)
2616  , data_values(data_values)
2617  {
2618  for (unsigned int d = 0; d < dim; ++d)
2619  {
2620  Assert(n_subintervals[d] >= 1,
2621  ExcMessage("There needs to be at least one subinterval in each "
2622  "coordinate direction."));
2623  Assert(interval_endpoints[d].first < interval_endpoints[d].second,
2624  ExcMessage("The interval in each coordinate direction needs "
2625  "to have positive size"));
2626  Assert(data_values.size()[d] == n_subintervals[d] + 1,
2627  ExcMessage("The data table does not have the correct size."));
2628  }
2629  }
2630 
2631 
2632  template <int dim>
2633  double
2634  InterpolatedUniformGridData<dim>::value(const Point<dim> & p,
2635  const unsigned int component) const
2636  {
2637  (void)component;
2638  Assert(
2639  component == 0,
2640  ExcMessage(
2641  "This is a scalar function object, the component can only be zero."));
2642 
2643  // find out where this data point lies, relative to the given
2644  // subdivision points
2645  TableIndices<dim> ix;
2646  for (unsigned int d = 0; d < dim; ++d)
2647  {
2648  const double delta_x =
2649  ((interval_endpoints[d].second - interval_endpoints[d].first) /
2650  n_subintervals[d]);
2651  if (p[d] <= interval_endpoints[d].first)
2652  ix[d] = 0;
2653  else if (p[d] >= interval_endpoints[d].second - delta_x)
2654  ix[d] = n_subintervals[d] - 1;
2655  else
2656  ix[d] =
2657  (unsigned int)((p[d] - interval_endpoints[d].first) / delta_x);
2658  }
2659 
2660  // now compute the relative point within the interval/rectangle/box
2661  // defined by the point coordinates found above. truncate below and
2662  // above to accommodate points that may lie outside the range
2663  Point<dim> p_unit;
2664  for (unsigned int d = 0; d < dim; ++d)
2665  {
2666  const double delta_x =
2667  ((interval_endpoints[d].second - interval_endpoints[d].first) /
2668  n_subintervals[d]);
2669  p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
2670  ix[d] * delta_x) /
2671  delta_x,
2672  1.),
2673  0.);
2674  }
2675 
2676  return interpolate(data_values, ix, p_unit);
2677  }
2678 
2679  /* ---------------------- Polynomial ----------------------- */
2680 
2681 
2682 
2683  template <int dim>
2684  Polynomial<dim>::Polynomial(const Table<2, double> & exponents,
2685  const std::vector<double> &coefficients)
2686  : Function<dim>(1)
2687  , exponents(exponents)
2688  , coefficients(coefficients)
2689  {
2690  Assert(exponents.n_rows() == coefficients.size(),
2691  ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
2692  Assert(exponents.n_cols() == dim,
2693  ExcDimensionMismatch(exponents.n_cols(), dim));
2694  }
2695 
2696 
2697 
2698  template <int dim>
2699  double
2700  Polynomial<dim>::value(const Point<dim> & p,
2701  const unsigned int component) const
2702  {
2703  (void)component;
2704  Assert(component == 0, ExcIndexRange(component, 0, 1));
2705 
2706  double sum = 0;
2707  for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2708  {
2709  double prod = 1;
2710  for (unsigned int s = 0; s < dim; ++s)
2711  {
2712  if (p[s] < 0)
2713  Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
2714  ExcMessage("Exponentiation of a negative base number with "
2715  "a real exponent can't be performed."));
2716  prod *= std::pow(p[s], exponents[monom][s]);
2717  }
2718  sum += coefficients[monom] * prod;
2719  }
2720  return sum;
2721  }
2722 
2723 
2724 
2725  template <int dim>
2726  void
2727  Polynomial<dim>::value_list(const std::vector<Point<dim>> &points,
2728  std::vector<double> & values,
2729  const unsigned int component) const
2730  {
2731  Assert(values.size() == points.size(),
2732  ExcDimensionMismatch(values.size(), points.size()));
2733 
2734  for (unsigned int i = 0; i < points.size(); ++i)
2735  values[i] = Polynomial<dim>::value(points[i], component);
2736  }
2737 
2738 
2739 
2740  template <int dim>
2741  Tensor<1, dim>
2742  Polynomial<dim>::gradient(const Point<dim> & p,
2743  const unsigned int component) const
2744  {
2745  (void)component;
2746  Assert(component == 0, ExcIndexRange(component, 0, 1));
2747 
2748  Tensor<1, dim> r;
2749 
2750  for (unsigned int d = 0; d < dim; ++d)
2751  {
2752  double sum = 0;
2753 
2754  for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
2755  {
2756  double prod = 1;
2757  for (unsigned int s = 0; s < dim; ++s)
2758  {
2759  if ((s == d) && (exponents[monom][s] == 0) && (p[s] == 0))
2760  {
2761  prod = 0;
2762  break;
2763  }
2764  else
2765  {
2766  if (p[s] < 0)
2767  Assert(std::floor(exponents[monom][s]) ==
2768  exponents[monom][s],
2769  ExcMessage(
2770  "Exponentiation of a negative base number with "
2771  "a real exponent can't be performed."));
2772  prod *=
2773  (s == d ? exponents[monom][s] *
2774  std::pow(p[s], exponents[monom][s] - 1) :
2775  std::pow(p[s], exponents[monom][s]));
2776  }
2777  }
2778  sum += coefficients[monom] * prod;
2779  }
2780  r[d] = sum;
2781  }
2782  return r;
2783  }
2784 
2785 
2786  // explicit instantiations
2787  template class SquareFunction<1>;
2788  template class SquareFunction<2>;
2789  template class SquareFunction<3>;
2790  template class Q1WedgeFunction<1>;
2791  template class Q1WedgeFunction<2>;
2792  template class Q1WedgeFunction<3>;
2793  template class PillowFunction<1>;
2794  template class PillowFunction<2>;
2795  template class PillowFunction<3>;
2796  template class CosineFunction<1>;
2797  template class CosineFunction<2>;
2798  template class CosineFunction<3>;
2799  template class CosineGradFunction<1>;
2800  template class CosineGradFunction<2>;
2801  template class CosineGradFunction<3>;
2802  template class ExpFunction<1>;
2803  template class ExpFunction<2>;
2804  template class ExpFunction<3>;
2805  template class JumpFunction<1>;
2806  template class JumpFunction<2>;
2807  template class JumpFunction<3>;
2808  template class FourierCosineFunction<1>;
2809  template class FourierCosineFunction<2>;
2810  template class FourierCosineFunction<3>;
2811  template class FourierSineFunction<1>;
2812  template class FourierSineFunction<2>;
2813  template class FourierSineFunction<3>;
2814  template class FourierCosineSum<1>;
2815  template class FourierCosineSum<2>;
2816  template class FourierCosineSum<3>;
2817  template class FourierSineSum<1>;
2818  template class FourierSineSum<2>;
2819  template class FourierSineSum<3>;
2820  template class SlitSingularityFunction<2>;
2821  template class SlitSingularityFunction<3>;
2822  template class Monomial<1>;
2823  template class Monomial<2>;
2824  template class Monomial<3>;
2825  template class Bessel1<1>;
2826  template class Bessel1<2>;
2827  template class Bessel1<3>;
2828  template class InterpolatedTensorProductGridData<1>;
2829  template class InterpolatedTensorProductGridData<2>;
2830  template class InterpolatedTensorProductGridData<3>;
2831  template class InterpolatedUniformGridData<1>;
2832  template class InterpolatedUniformGridData<2>;
2833  template class InterpolatedUniformGridData<3>;
2834  template class Polynomial<1>;
2835  template class Polynomial<2>;
2836  template class Polynomial<3>;
2837 } // namespace Functions
2838 
2839 DEAL_II_NAMESPACE_CLOSE
Functions::Monomial::exponents
const Tensor< 1, dim > exponents
Definition: function_lib.h:1185
Functions::FourierCosineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2136
Functions::CosineGradFunction::value
virtual double value(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:726
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:1366
Function< dim >::n_components
const unsigned int n_components
Definition: function.h:160
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:654
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:786
Functions::Monomial::Monomial
Monomial(const Tensor< 1, dim > &exponents, const unsigned int n_components=1)
Definition: function_lib.cc:2176
Functions::Polynomial::Polynomial
Polynomial(const Table< 2, double > &exponents, const std::vector< double > &coefficients)
Definition: function_lib.cc:2684
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
Vector::size
size_type size() const
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1407
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
AssertVectorVectorDimension
#define AssertVectorVectorDimension(vec, dim1, dim2)
Definition: exceptions.h:1377
Functions::Monomial::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:2208
Functions::FourierSineSum::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2062
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Functions::FourierCosineSum::FourierCosineSum
FourierCosineSum(const std::vector< Point< dim >> &fourier_coefficients, const std::vector< double > &weights)
Definition: function_lib.cc:2102
Functions::FourierCosineFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1962
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:1988
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:2186
Functions::JumpFunction::memory_consumption
std::size_t memory_consumption() const
Definition: function_lib.cc:1915
Functions::InterpolatedTensorProductGridData::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2550
Functions::FourierSineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2000
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:1072
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:1938
Functions::Monomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2221
Functions::FourierCosineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2154
numbers::PI
static const double PI
Definition: numbers.h:143
Functions::InterpolatedTensorProductGridData::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2581
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:2285
Functions::SlitSingularityFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1557
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:2118
Functions::FourierCosineFunction::fourier_coefficients
const Tensor< 1, dim > fourier_coefficients
Definition: function_lib.h:743
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:2260
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:1119
Functions::JumpFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1881
TableBase::size
size_type size(const unsigned int i) const
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:1227
Functions::FourierSineSum::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2080
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:1533
Functions::InterpolatedUniformGridData::data_values
const Table< dim, double > data_values
Definition: function_lib.h:1375
Functions::Bessel1::Bessel1
Bessel1(const unsigned int order, const double wave_number, const Point< dim > center=Point< dim >())
Definition: function_lib.cc:2273
Functions::CosineGradFunction::vector_value
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
Definition: function_lib.cc:752
Functions::SlitSingularityFunction::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1470
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:905
Functions::InterpolatedTensorProductGridData::table_index_of_point
TableIndices< dim > table_index_of_point(const Point< dim > &p) const
Definition: function_lib.cc:2513
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:1144
Functions::CosineGradFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:872
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:1043
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:1433
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:1860
Functions::Polynomial::coefficients
const std::vector< double > coefficients
Definition: function_lib.h:1438
Point::square
numbers::NumberTraits< Number >::real_type square() const
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:822
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:1508
Functions::JumpFunction::direction
const Point< dim > direction
Definition: function_lib.h:668
Functions::FourierSineSum::fourier_coefficients
const std::vector< Point< dim > > fourier_coefficients
Definition: function_lib.h:850
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:2295
Functions::FourierCosineFunction::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1950
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:1823
Functions::PillowFunction::PillowFunction
PillowFunction(const double offset=0.)
Definition: function_lib.cc:254
Functions::InterpolatedTensorProductGridData::coordinate_values
const std::array< std::vector< double >, dim > coordinate_values
Definition: function_lib.h:1281
numbers::PI_2
static const double PI_2
Definition: numbers.h:148
Functions::InterpolatedTensorProductGridData::InterpolatedTensorProductGridData
InterpolatedTensorProductGridData(const std::array< std::vector< double >, dim > &coordinate_values, const Table< dim, double > &data_values)
Definition: function_lib.cc:2485
Functions::InterpolatedUniformGridData::n_subintervals
const std::array< unsigned int, dim > n_subintervals
Definition: function_lib.h:1370
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:2610
Functions::FourierCosineFunction::FourierCosineFunction
FourierCosineFunction(const Tensor< 1, dim > &fourier_coefficients)
Definition: function_lib.cc:1928
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:1025
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:2028
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:1978
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:2634
Functions::FourierSineFunction::fourier_coefficients
const Tensor< 1, dim > fourier_coefficients
Definition: function_lib.h:798
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:2700
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:1485
Functions::Polynomial::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2742
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:1895
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:1090
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:1365
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:1833
Functions::CosineGradFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component) const override
Definition: function_lib.cc:863
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:1542
Functions::FourierSineSum::value
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2044
Functions::JumpFunction::laplacian
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:1850
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:2012
Functions::SquareFunction::laplacian_list
virtual void laplacian_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
Definition: function_lib.cc:77
Functions::Bessel1::gradient
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: function_lib.cc:2311
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:2727
Functions::JumpFunction::JumpFunction
JumpFunction(const Point< dim > &direction, const double steepness)
Definition: function_lib.cc:1799
Generated by   doxygen 1.8.11