Reference documentation for deal.II version Git e159ac89de 2020-06-06 19:38:41 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
function_lib_cutoff.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2020 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 
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 
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  {
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  {
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  {
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  {
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  {
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  AssertIndexRange(component, this->n_components);
255 
256 
258  component == this->selected)
259  for (unsigned int k = 0; k < values.size(); ++k)
260  values[k] = (this->center.distance(points[k]) < this->radius) ?
261  this->rescaling :
262  0.;
263  else
264  std::fill(values.begin(), values.end(), 0.);
265  }
266 
267 
268  template <int dim>
269  void
271  const std::vector<Point<dim>> &points,
272  std::vector<Vector<double>> & values) const
273  {
274  Assert(values.size() == points.size(),
275  ExcDimensionMismatch(values.size(), points.size()));
276 
277  for (unsigned int k = 0; k < values.size(); ++k)
278  {
279  const double val = (this->center.distance(points[k]) < this->radius) ?
280  this->rescaling :
281  0.;
283  values[k] = val;
284  else
285  {
286  values[k] = 0;
287  values[k](this->selected) = val;
288  }
289  }
290  }
291 
292  template <int dim>
294  const Point<dim> p,
295  const unsigned int n_components,
296  const unsigned int select,
297  const bool integrate_to_one)
298  : CutOffFunctionBase<dim>(r,
299  p,
300  n_components,
301  select,
302  integrate_to_one,
303  integral_W1[dim - 1])
304  {}
305 
306 
307  template <int dim>
308  double
310  const unsigned int component) const
311  {
313  component == this->selected)
314  {
315  const double d = this->center.distance(p);
316  return ((d < this->radius) ?
317  (this->radius - d) / this->radius * this->rescaling :
318  0.);
319  }
320  return 0.;
321  }
322 
323 
324  template <int dim>
325  void
326  CutOffFunctionW1<dim>::value_list(const std::vector<Point<dim>> &points,
327  std::vector<double> & values,
328  const unsigned int component) const
329  {
330  Assert(values.size() == points.size(),
331  ExcDimensionMismatch(values.size(), points.size()));
332 
334  component == this->selected)
335  for (unsigned int i = 0; i < values.size(); ++i)
336  {
337  const double d = this->center.distance(points[i]);
338  values[i] = ((d < this->radius) ?
339  (this->radius - d) / this->radius * this->rescaling :
340  0.);
341  }
342  else
343  std::fill(values.begin(), values.end(), 0.);
344  }
345 
346 
347 
348  template <int dim>
349  void
351  const std::vector<Point<dim>> &points,
352  std::vector<Vector<double>> & values) const
353  {
354  Assert(values.size() == points.size(),
355  ExcDimensionMismatch(values.size(), points.size()));
356 
357  for (unsigned int k = 0; k < values.size(); ++k)
358  {
359  const double d = this->center.distance(points[k]);
360  const double val =
361  (d < this->radius) ?
362  (this->radius - d) / this->radius * this->rescaling :
363  0.;
365  values[k] = val;
366  else
367  {
368  values[k] = 0;
369  values[k](this->selected) = val;
370  }
371  }
372  }
373 
374 
375  template <int dim>
377  const double r,
378  const Point<dim> p,
379  const unsigned int n_components,
380  const unsigned int select,
381  bool integrate_to_one)
382  : CutOffFunctionBase<dim>(r,
383  p,
384  n_components,
385  select,
386  integrate_to_one,
387  integral_Cinfty[dim - 1])
388  {}
389 
390 
391  template <int dim>
392  double
394  const unsigned int component) const
395  {
397  component == this->selected)
398  {
399  const double d = this->center.distance(p);
400  const double r = this->radius;
401  if (d >= r)
402  return 0.;
403  const double e = -r * r / (r * r - d * d);
404  return ((e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling);
405  }
406  return 0.;
407  }
408 
409 
410  template <int dim>
411  void
413  std::vector<double> & values,
414  const unsigned int component) const
415  {
416  Assert(values.size() == points.size(),
417  ExcDimensionMismatch(values.size(), points.size()));
418 
419  const double r = this->radius;
420 
422  component == this->selected)
423  for (unsigned int i = 0; i < values.size(); ++i)
424  {
425  const double d = this->center.distance(points[i]);
426  if (d >= r)
427  {
428  values[i] = 0.;
429  }
430  else
431  {
432  const double e = -r * r / (r * r - d * d);
433  values[i] =
434  (e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling;
435  }
436  }
437  else
438  std::fill(values.begin(), values.end(), 0.);
439  }
440 
441 
442  template <int dim>
443  void
445  const std::vector<Point<dim>> &points,
446  std::vector<Vector<double>> & values) const
447  {
448  Assert(values.size() == points.size(),
449  ExcDimensionMismatch(values.size(), points.size()));
450 
451  for (unsigned int k = 0; k < values.size(); ++k)
452  {
453  const double d = this->center.distance(points[k]);
454  const double r = this->radius;
455  double val = 0.;
456  if (d < this->radius)
457  {
458  const double e = -r * r / (r * r - d * d);
459  if (e > -50)
460  val = numbers::E * std::exp(e) * this->rescaling;
461  }
462 
464  values[k] = val;
465  else
466  {
467  values[k] = 0;
468  values[k](this->selected) = val;
469  }
470  }
471  }
472 
473 
474 
475  template <int dim>
478  const unsigned int) const
479  {
480  const double d = this->center.distance(p);
481  const double r = this->radius;
482  if (d >= r)
483  return Tensor<1, dim>();
484  const double e = -d * d / (r - d) / (r + d);
485  return ((e < -50) ? Point<dim>() :
486  (p - this->center) / d *
487  (-2.0 * r * r / std::pow(-r * r + d * d, 2.0) * d *
488  std::exp(e)) *
489  this->rescaling);
490  }
491 
492 
493 
494  template <int dim>
496  const Point<dim> p,
497  const unsigned int n_components,
498  const unsigned int select,
499  bool integrate_to_one)
500  : CutOffFunctionBase<dim>(r,
501  p,
502  n_components,
503  select,
504  integrate_to_one,
505  integral_C1[dim - 1])
506  {}
507 
508 
509  template <int dim>
510  double
512  const unsigned int component) const
513  {
515  component == this->selected)
516  {
517  const double d = this->center.distance(p);
518  const double r = this->radius;
519  if (d >= r)
520  return 0.;
521  return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
522  }
523  return 0.;
524  }
525 
526 
527  template <int dim>
528  void
529  CutOffFunctionC1<dim>::value_list(const std::vector<Point<dim>> &points,
530  std::vector<double> & values,
531  const unsigned int component) const
532  {
533  Assert(values.size() == points.size(),
534  ExcDimensionMismatch(values.size(), points.size()));
535 
536  const double r = this->radius;
537 
539  component == this->selected)
540  for (unsigned int i = 0; i < values.size(); ++i)
541  {
542  const double d = this->center.distance(points[i]);
543  if (d >= r)
544  {
545  values[i] = 0.;
546  }
547  else
548  {
549  values[i] =
550  .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
551  }
552  }
553  else
554  std::fill(values.begin(), values.end(), 0.);
555  }
556 
557 
558  template <int dim>
559  void
561  const std::vector<Point<dim>> &points,
562  std::vector<Vector<double>> & values) const
563  {
564  Assert(values.size() == points.size(),
565  ExcDimensionMismatch(values.size(), points.size()));
566 
567  for (unsigned int k = 0; k < values.size(); ++k)
568  {
569  const double d = this->center.distance(points[k]);
570  const double r = this->radius;
571  double val = 0.;
572  if (d < this->radius)
573  {
574  val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
575  }
576 
578  values[k] = val;
579  else
580  {
581  values[k] = 0;
582  values[k](this->selected) = val;
583  }
584  }
585  }
586 
587 
588 
589  template <int dim>
591  CutOffFunctionC1<dim>::gradient(const Point<dim> &p, const unsigned int) const
592  {
593  const double d = this->center.distance(p);
594  const double r = this->radius;
595  if (d >= r)
596  return Tensor<1, dim>();
597  return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) *
598  (p - this->center) / d * this->rescaling;
599  }
600 
601 
602  // explicit instantiations
603  template class CutOffFunctionBase<1>;
604  template class CutOffFunctionBase<2>;
605  template class CutOffFunctionBase<3>;
606 
607  template class CutOffFunctionLinfty<1>;
608  template class CutOffFunctionLinfty<2>;
609  template class CutOffFunctionLinfty<3>;
610 
611  template class CutOffFunctionW1<1>;
612  template class CutOffFunctionW1<2>;
613  template class CutOffFunctionW1<3>;
614 
615  template class CutOffFunctionCinfty<1>;
616  template class CutOffFunctionCinfty<2>;
617  template class CutOffFunctionCinfty<3>;
618 
619  template class CutOffFunctionC1<1>;
620  template class CutOffFunctionC1<2>;
621  template class CutOffFunctionC1<3>;
622 
623  template class CutOffFunctionTensorProduct<1>;
624  template class CutOffFunctionTensorProduct<2>;
625  template class CutOffFunctionTensorProduct<3>;
626 } // namespace Functions
627 
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
const unsigned int n_components
Definition: function.h:164
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
virtual void set_center(const Point< dim > &center) override
static ::ExceptionBase & ExcNotInitialized()
std::array< std::unique_ptr< CutOffFunctionBase< 1 > >, dim > base
virtual void set_center(const Point< dim > &p)
T fixed_power(const T t)
Definition: utilities.h:1045
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:206
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void set_radius(const double radius) override
#define Assert(cond, exc)
Definition: exceptions.h:1403
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
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
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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
const unsigned int selected
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:231
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
Point< 3 > center
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)