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