Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
function_lib_cutoff.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/function_lib.h>
17 #include <deal.II/base/point.h>
18 #include <deal.II/base/tensor.h>
19 
20 #include <deal.II/lac/vector.h>
21 
22 #include <cmath>
23 
24 DEAL_II_NAMESPACE_OPEN
25 
26 
27 namespace Functions
28 {
29  template <int dim>
31  const double r,
32  const Point<dim> p,
33  const unsigned int n_components,
34  const unsigned int select,
35  const bool integrate_to_one,
36  const double unitary_integral_value)
37  : Function<dim>(n_components)
38  , center(p)
39  , radius(r)
40  , selected(select)
41  , integrate_to_one(integrate_to_one)
42  , unitary_integral_value(unitary_integral_value)
43  , rescaling(integrate_to_one ? 1. / (unitary_integral_value *
44  Utilities::fixed_power<dim>(radius)) :
45  1.0)
46  {
47  Assert(r > 0, ExcMessage("You must specify a radius > 0."));
48  }
49 
50 
51  template <int dim>
52  void
54  {
55  set_center(p);
56  }
57 
58 
59 
60  template <int dim>
61  void
63  {
64  center = p;
65  }
66 
67 
68 
69  template <int dim>
70  const Point<dim> &
72  {
73  return center;
74  }
75 
76 
77 
78  template <int dim>
79  void
81  {
82  set_radius(r);
83  }
84 
85 
86 
87  template <int dim>
88  void
90  {
91  radius = r;
92  Assert(r > 0, ExcMessage("You must specify a radius > 0."));
93  if (integrate_to_one)
94  rescaling =
95  1. / (unitary_integral_value * Utilities::fixed_power<dim>(radius));
96  else
97  rescaling = 1.0;
98  }
99 
100 
101 
102  template <int dim>
103  double
105  {
106  return radius;
107  }
108 
109 
110 
111  template <int dim>
112  bool
114  {
115  return integrate_to_one;
116  }
117 
118 
119 
120  template <int dim>
122  double radius,
123  const Point<dim> & center,
124  const unsigned int n_components,
125  const unsigned int select,
126  const bool integrate_to_one)
127  : CutOffFunctionBase<dim>(radius,
128  center,
129  n_components,
130  select,
131  integrate_to_one)
132  , initialized(false)
133  {}
134 
135 
136 
137  template <int dim>
138  void
140  {
141  Assert(initialized, ExcNotInitialized());
142  for (unsigned int i = 0; i < dim; ++i)
143  base[i]->set_center(Point<1>(p[i]));
145  }
146 
147 
148 
149  template <int dim>
150  void
152  {
153  Assert(initialized, ExcNotInitialized());
154  for (unsigned int i = 0; i < dim; ++i)
155  base[i]->set_radius(r);
157  }
158 
159 
160 
161  template <int dim>
162  double
164  const unsigned int component) const
165  {
166  Assert(initialized, ExcNotInitialized());
167  double ret = 1.0;
168  for (unsigned int i = 0; i < dim; ++i)
169  ret *= base[i]->value(Point<1>(p[i]), component);
170  return ret;
171  }
172 
173 
174 
175  template <int dim>
178  const unsigned int component) const
179  {
180  Assert(initialized, ExcNotInitialized());
181  Tensor<1, dim> ret;
182  for (unsigned int i = 0; i < dim; ++i)
183  ret[i] = base[i]->gradient(Point<1>(p[i]), component)[0];
184  return ret;
185  }
186 
187 
188 
190  namespace
191  {
192  // Integral of CutOffFunctionLinfty in dimension 1, 2, and 3 when the radius
193  // is one
194  const double integral_Linfty[] = {2.0,
195  3.14159265358979323846264338328,
196  4.18879020478639098461685784437};
197 
198  // Integral of CutOffFunctionW1 in dimension 1, 2, and 3 when the radius
199  // is one
200  const double integral_W1[] = {1.0,
201  1.04719755119659774615421446109,
202  1.04719755119659774615421446109};
203 
204  // Integral of CutOffFunctionCinfty in dimension 1, 2, and 3 when the radius
205  // is one
206  const double integral_Cinfty[] = {1.20690032243787617533623799633,
207  1.26811216112759608094632335664,
208  1.1990039070192139033798473858};
209 
210  // Integral of CutOffFunctionC1 in dimension 1, 2, and 3 when the radius
211  // is one
212  const double integral_C1[] = {1.0,
213  0.93417655442731527615578663815,
214  0.821155557658032806157358815206};
215  } // namespace
216 
217 
218  template <int dim>
220  const double r,
221  const Point<dim> p,
222  const unsigned int n_components,
223  const unsigned int select,
224  const bool integrate_to_one)
225  : CutOffFunctionBase<dim>(r,
226  p,
227  n_components,
228  select,
229  integrate_to_one,
230  integral_Linfty[dim - 1])
231  {}
232 
233 
234  template <int dim>
235  double
237  const unsigned int component) const
238  {
239  if (this->selected == CutOffFunctionBase<dim>::no_component ||
240  component == this->selected)
241  return ((this->center.distance(p) < this->radius) ? this->rescaling : 0.);
242  return 0.;
243  }
244 
245 
246  template <int dim>
247  void
249  std::vector<double> & values,
250  const unsigned int component) const
251  {
252  Assert(values.size() == points.size(),
253  ExcDimensionMismatch(values.size(), points.size()));
254  Assert(component < this->n_components,
255  ExcIndexRange(component, 0, this->n_components));
256 
257 
258  if (this->selected == CutOffFunctionBase<dim>::no_component ||
259  component == this->selected)
260  for (unsigned int k = 0; k < values.size(); ++k)
261  values[k] = (this->center.distance(points[k]) < this->radius) ?
262  this->rescaling :
263  0.;
264  else
265  std::fill(values.begin(), values.end(), 0.);
266  }
267 
268 
269  template <int dim>
270  void
272  const std::vector<Point<dim>> &points,
273  std::vector<Vector<double>> & values) const
274  {
275  Assert(values.size() == points.size(),
276  ExcDimensionMismatch(values.size(), points.size()));
277 
278  for (unsigned int k = 0; k < values.size(); ++k)
279  {
280  const double val = (this->center.distance(points[k]) < this->radius) ?
281  this->rescaling :
282  0.;
283  if (this->selected == CutOffFunctionBase<dim>::no_component)
284  values[k] = val;
285  else
286  {
287  values[k] = 0;
288  values[k](this->selected) = val;
289  }
290  }
291  }
292 
293  template <int dim>
295  const Point<dim> p,
296  const unsigned int n_components,
297  const unsigned int select,
298  const bool integrate_to_one)
299  : CutOffFunctionBase<dim>(r,
300  p,
301  n_components,
302  select,
303  integrate_to_one,
304  integral_W1[dim - 1])
305  {}
306 
307 
308  template <int dim>
309  double
311  const unsigned int component) const
312  {
313  if (this->selected == CutOffFunctionBase<dim>::no_component ||
314  component == this->selected)
315  {
316  const double d = this->center.distance(p);
317  return ((d < this->radius) ?
318  (this->radius - d) / this->radius * this->rescaling :
319  0.);
320  }
321  return 0.;
322  }
323 
324 
325  template <int dim>
326  void
327  CutOffFunctionW1<dim>::value_list(const std::vector<Point<dim>> &points,
328  std::vector<double> & values,
329  const unsigned int component) const
330  {
331  Assert(values.size() == points.size(),
332  ExcDimensionMismatch(values.size(), points.size()));
333 
334  if (this->selected == CutOffFunctionBase<dim>::no_component ||
335  component == this->selected)
336  for (unsigned int i = 0; i < values.size(); ++i)
337  {
338  const double d = this->center.distance(points[i]);
339  values[i] = ((d < this->radius) ?
340  (this->radius - d) / this->radius * this->rescaling :
341  0.);
342  }
343  else
344  std::fill(values.begin(), values.end(), 0.);
345  }
346 
347 
348 
349  template <int dim>
350  void
352  const std::vector<Point<dim>> &points,
353  std::vector<Vector<double>> & values) const
354  {
355  Assert(values.size() == points.size(),
356  ExcDimensionMismatch(values.size(), points.size()));
357 
358  for (unsigned int k = 0; k < values.size(); ++k)
359  {
360  const double d = this->center.distance(points[k]);
361  const double val =
362  (d < this->radius) ?
363  (this->radius - d) / this->radius * this->rescaling :
364  0.;
365  if (this->selected == CutOffFunctionBase<dim>::no_component)
366  values[k] = val;
367  else
368  {
369  values[k] = 0;
370  values[k](this->selected) = val;
371  }
372  }
373  }
374 
375 
376  template <int dim>
378  const double r,
379  const Point<dim> p,
380  const unsigned int n_components,
381  const unsigned int select,
382  bool integrate_to_one)
383  : CutOffFunctionBase<dim>(r,
384  p,
385  n_components,
386  select,
387  integrate_to_one,
388  integral_Cinfty[dim - 1])
389  {}
390 
391 
392  template <int dim>
393  double
395  const unsigned int component) const
396  {
397  if (this->selected == CutOffFunctionBase<dim>::no_component ||
398  component == this->selected)
399  {
400  const double d = this->center.distance(p);
401  const double r = this->radius;
402  if (d >= r)
403  return 0.;
404  const double e = -r * r / (r * r - d * d);
405  return ((e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling);
406  }
407  return 0.;
408  }
409 
410 
411  template <int dim>
412  void
414  std::vector<double> & values,
415  const unsigned int component) const
416  {
417  Assert(values.size() == points.size(),
418  ExcDimensionMismatch(values.size(), points.size()));
419 
420  const double r = this->radius;
421 
422  if (this->selected == CutOffFunctionBase<dim>::no_component ||
423  component == this->selected)
424  for (unsigned int i = 0; i < values.size(); ++i)
425  {
426  const double d = this->center.distance(points[i]);
427  if (d >= r)
428  {
429  values[i] = 0.;
430  }
431  else
432  {
433  const double e = -r * r / (r * r - d * d);
434  values[i] =
435  (e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling;
436  }
437  }
438  else
439  std::fill(values.begin(), values.end(), 0.);
440  }
441 
442 
443  template <int dim>
444  void
446  const std::vector<Point<dim>> &points,
447  std::vector<Vector<double>> & values) const
448  {
449  Assert(values.size() == points.size(),
450  ExcDimensionMismatch(values.size(), points.size()));
451 
452  for (unsigned int k = 0; k < values.size(); ++k)
453  {
454  const double d = this->center.distance(points[k]);
455  const double r = this->radius;
456  double val = 0.;
457  if (d < this->radius)
458  {
459  const double e = -r * r / (r * r - d * d);
460  if (e > -50)
461  val = numbers::E * std::exp(e) * this->rescaling;
462  }
463 
464  if (this->selected == CutOffFunctionBase<dim>::no_component)
465  values[k] = val;
466  else
467  {
468  values[k] = 0;
469  values[k](this->selected) = val;
470  }
471  }
472  }
473 
474 
475 
476  template <int dim>
479  const unsigned int) const
480  {
481  const double d = this->center.distance(p);
482  const double r = this->radius;
483  if (d >= r)
484  return Tensor<1, dim>();
485  const double e = -d * d / (r - d) / (r + d);
486  return ((e < -50) ? Point<dim>() :
487  (p - this->center) / d *
488  (-2.0 * r * r / std::pow(-r * r + d * d, 2.0) * d *
489  std::exp(e)) *
490  this->rescaling);
491  }
492 
493 
494 
495  template <int dim>
497  const Point<dim> p,
498  const unsigned int n_components,
499  const unsigned int select,
500  bool integrate_to_one)
501  : CutOffFunctionBase<dim>(r,
502  p,
503  n_components,
504  select,
505  integrate_to_one,
506  integral_C1[dim - 1])
507  {}
508 
509 
510  template <int dim>
511  double
513  const unsigned int component) const
514  {
515  if (this->selected == CutOffFunctionBase<dim>::no_component ||
516  component == this->selected)
517  {
518  const double d = this->center.distance(p);
519  const double r = this->radius;
520  if (d >= r)
521  return 0.;
522  return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
523  }
524  return 0.;
525  }
526 
527 
528  template <int dim>
529  void
530  CutOffFunctionC1<dim>::value_list(const std::vector<Point<dim>> &points,
531  std::vector<double> & values,
532  const unsigned int component) const
533  {
534  Assert(values.size() == points.size(),
535  ExcDimensionMismatch(values.size(), points.size()));
536 
537  const double r = this->radius;
538 
539  if (this->selected == CutOffFunctionBase<dim>::no_component ||
540  component == this->selected)
541  for (unsigned int i = 0; i < values.size(); ++i)
542  {
543  const double d = this->center.distance(points[i]);
544  if (d >= r)
545  {
546  values[i] = 0.;
547  }
548  else
549  {
550  values[i] =
551  .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
552  }
553  }
554  else
555  std::fill(values.begin(), values.end(), 0.);
556  }
557 
558 
559  template <int dim>
560  void
562  const std::vector<Point<dim>> &points,
563  std::vector<Vector<double>> & values) const
564  {
565  Assert(values.size() == points.size(),
566  ExcDimensionMismatch(values.size(), points.size()));
567 
568  for (unsigned int k = 0; k < values.size(); ++k)
569  {
570  const double d = this->center.distance(points[k]);
571  const double r = this->radius;
572  double val = 0.;
573  if (d < this->radius)
574  {
575  val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
576  }
577 
578  if (this->selected == CutOffFunctionBase<dim>::no_component)
579  values[k] = val;
580  else
581  {
582  values[k] = 0;
583  values[k](this->selected) = val;
584  }
585  }
586  }
587 
588 
589 
590  template <int dim>
592  CutOffFunctionC1<dim>::gradient(const Point<dim> &p, const unsigned int) const
593  {
594  const double d = this->center.distance(p);
595  const double r = this->radius;
596  if (d >= r)
597  return Tensor<1, dim>();
598  return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) *
599  (p - this->center) / d * this->rescaling;
600  }
601 
602 
603  // explicit instantiations
604  template class CutOffFunctionBase<1>;
605  template class CutOffFunctionBase<2>;
606  template class CutOffFunctionBase<3>;
607 
608  template class CutOffFunctionLinfty<1>;
609  template class CutOffFunctionLinfty<2>;
610  template class CutOffFunctionLinfty<3>;
611 
612  template class CutOffFunctionW1<1>;
613  template class CutOffFunctionW1<2>;
614  template class CutOffFunctionW1<3>;
615 
616  template class CutOffFunctionCinfty<1>;
617  template class CutOffFunctionCinfty<2>;
618  template class CutOffFunctionCinfty<3>;
619 
620  template class CutOffFunctionC1<1>;
621  template class CutOffFunctionC1<2>;
622  template class CutOffFunctionC1<3>;
623 
624  template class CutOffFunctionTensorProduct<1>;
625  template class CutOffFunctionTensorProduct<2>;
626  template class CutOffFunctionTensorProduct<3>;
627 } // namespace Functions
628 
629 DEAL_II_NAMESPACE_CLOSE
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
CutOffFunctionW1(const double radius=1., const Point< dim >=Point< dim >(), const unsigned int n_components=1, const unsigned int select=CutOffFunctionBase< dim >::no_component, const bool integrate_to_one=false)
CutOffFunctionC1(const double radius=1., const Point< dim >=Point< dim >(), const unsigned int n_components=1, const unsigned int select=CutOffFunctionBase< dim >::no_component, bool integrate_to_one=false)
virtual void set_center(const Point< dim > &center) override
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void set_center(const Point< dim > &p)
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
CutOffFunctionBase(const double radius=1., const Point< dim > center=Point< dim >(), const unsigned int n_components=1, const unsigned int select=CutOffFunctionBase< dim >::no_component, const bool integrate_to_one=false, const double unitary_integral_value=1.0)
static constexpr double E
Definition: numbers.h:121
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void set_radius(const double radius) override
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
CutOffFunctionLinfty(const double radius=1., const Point< dim >=Point< dim >(), const unsigned int n_components=1, const unsigned int select=CutOffFunctionBase< dim >::no_component, const bool integrate_to_one=false)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
Definition: cuda.h:31
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
static constexpr double PI
Definition: numbers.h:146
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
CutOffFunctionCinfty(const double radius=1., const Point< dim >=Point< dim >(), const unsigned int n_components=1, const unsigned int select=CutOffFunctionBase< dim >::no_component, bool integrate_to_one=false)
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual void set_radius(const double r)