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