Reference documentation for deal.II version Git ede8f93e86 2020-12-03 14:59:20 -0700
\(\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 
52  template <int dim>
53  void
55  {
56  center = p;
57  }
58 
59 
60 
61  template <int dim>
62  const Point<dim> &
64  {
65  return center;
66  }
67 
68 
69 
70  template <int dim>
71  void
73  {
74  radius = r;
75  Assert(r > 0, ExcMessage("You must specify a radius > 0."));
76  if (integrate_to_one)
77  rescaling =
78  1. / (unitary_integral_value * Utilities::fixed_power<dim>(radius));
79  else
80  rescaling = 1.0;
81  }
82 
83 
84 
85  template <int dim>
86  double
88  {
89  return radius;
90  }
91 
92 
93 
94  template <int dim>
95  bool
97  {
98  return integrate_to_one;
99  }
100 
101 
102 
103  template <int dim>
105  double radius,
106  const Point<dim> & center,
107  const unsigned int n_components,
108  const unsigned int select,
109  const bool integrate_to_one)
110  : CutOffFunctionBase<dim>(radius,
111  center,
112  n_components,
113  select,
114  integrate_to_one)
115  , initialized(false)
116  {}
117 
118 
119 
120  template <int dim>
121  void
123  {
125  for (unsigned int i = 0; i < dim; ++i)
126  base[i]->set_center(Point<1>(p[i]));
128  }
129 
130 
131 
132  template <int dim>
133  void
135  {
137  for (unsigned int i = 0; i < dim; ++i)
138  base[i]->set_radius(r);
140  }
141 
142 
143 
144  template <int dim>
145  double
147  const unsigned int component) const
148  {
150  double ret = 1.0;
151  for (unsigned int i = 0; i < dim; ++i)
152  ret *= base[i]->value(Point<1>(p[i]), component);
153  return ret;
154  }
155 
156 
157 
158  template <int dim>
161  const unsigned int component) const
162  {
164  Tensor<1, dim> ret;
165  for (unsigned int i = 0; i < dim; ++i)
166  ret[i] = base[i]->gradient(Point<1>(p[i]), component)[0];
167  return ret;
168  }
169 
170 
171 
173  namespace
174  {
175  // Integral of CutOffFunctionLinfty in dimension 1, 2, and 3 when the radius
176  // is one
177  const double integral_Linfty[] = {2.0,
178  3.14159265358979323846264338328,
179  4.18879020478639098461685784437};
180 
181  // Integral of CutOffFunctionW1 in dimension 1, 2, and 3 when the radius
182  // is one
183  const double integral_W1[] = {1.0,
184  1.04719755119659774615421446109,
185  1.04719755119659774615421446109};
186 
187  // Integral of CutOffFunctionCinfty in dimension 1, 2, and 3 when the radius
188  // is one
189  const double integral_Cinfty[] = {1.20690032243787617533623799633,
190  1.26811216112759608094632335664,
191  1.1990039070192139033798473858};
192 
193  // Integral of CutOffFunctionC1 in dimension 1, 2, and 3 when the radius
194  // is one
195  const double integral_C1[] = {1.0,
196  0.93417655442731527615578663815,
197  0.821155557658032806157358815206};
198  } // namespace
199 
200 
201  template <int dim>
203  const double r,
204  const Point<dim> p,
205  const unsigned int n_components,
206  const unsigned int select,
207  const bool integrate_to_one)
208  : CutOffFunctionBase<dim>(r,
209  p,
210  n_components,
211  select,
212  integrate_to_one,
213  integral_Linfty[dim - 1])
214  {}
215 
216 
217  template <int dim>
218  double
220  const unsigned int component) const
221  {
223  component == this->selected)
224  return ((this->center.distance(p) < this->radius) ? this->rescaling : 0.);
225  return 0.;
226  }
227 
228 
229  template <int dim>
230  void
232  std::vector<double> & values,
233  const unsigned int component) const
234  {
235  Assert(values.size() == points.size(),
236  ExcDimensionMismatch(values.size(), points.size()));
237  AssertIndexRange(component, this->n_components);
238 
239 
241  component == this->selected)
242  for (unsigned int k = 0; k < values.size(); ++k)
243  values[k] = (this->center.distance(points[k]) < this->radius) ?
244  this->rescaling :
245  0.;
246  else
247  std::fill(values.begin(), values.end(), 0.);
248  }
249 
250 
251  template <int dim>
252  void
254  const std::vector<Point<dim>> &points,
255  std::vector<Vector<double>> & values) const
256  {
257  Assert(values.size() == points.size(),
258  ExcDimensionMismatch(values.size(), points.size()));
259 
260  for (unsigned int k = 0; k < values.size(); ++k)
261  {
262  const double val = (this->center.distance(points[k]) < this->radius) ?
263  this->rescaling :
264  0.;
266  values[k] = val;
267  else
268  {
269  values[k] = 0;
270  values[k](this->selected) = val;
271  }
272  }
273  }
274 
275  template <int dim>
277  const Point<dim> p,
278  const unsigned int n_components,
279  const unsigned int select,
280  const bool integrate_to_one)
281  : CutOffFunctionBase<dim>(r,
282  p,
283  n_components,
284  select,
285  integrate_to_one,
286  integral_W1[dim - 1])
287  {}
288 
289 
290  template <int dim>
291  double
293  const unsigned int component) const
294  {
296  component == this->selected)
297  {
298  const double d = this->center.distance(p);
299  return ((d < this->radius) ?
300  (this->radius - d) / this->radius * this->rescaling :
301  0.);
302  }
303  return 0.;
304  }
305 
306 
307  template <int dim>
308  void
309  CutOffFunctionW1<dim>::value_list(const std::vector<Point<dim>> &points,
310  std::vector<double> & values,
311  const unsigned int component) const
312  {
313  Assert(values.size() == points.size(),
314  ExcDimensionMismatch(values.size(), points.size()));
315 
317  component == this->selected)
318  for (unsigned int i = 0; i < values.size(); ++i)
319  {
320  const double d = this->center.distance(points[i]);
321  values[i] = ((d < this->radius) ?
322  (this->radius - d) / this->radius * this->rescaling :
323  0.);
324  }
325  else
326  std::fill(values.begin(), values.end(), 0.);
327  }
328 
329 
330 
331  template <int dim>
332  void
334  const std::vector<Point<dim>> &points,
335  std::vector<Vector<double>> & values) const
336  {
337  Assert(values.size() == points.size(),
338  ExcDimensionMismatch(values.size(), points.size()));
339 
340  for (unsigned int k = 0; k < values.size(); ++k)
341  {
342  const double d = this->center.distance(points[k]);
343  const double val =
344  (d < this->radius) ?
345  (this->radius - d) / this->radius * this->rescaling :
346  0.;
348  values[k] = val;
349  else
350  {
351  values[k] = 0;
352  values[k](this->selected) = val;
353  }
354  }
355  }
356 
357 
358  template <int dim>
360  const double r,
361  const Point<dim> p,
362  const unsigned int n_components,
363  const unsigned int select,
364  bool integrate_to_one)
365  : CutOffFunctionBase<dim>(r,
366  p,
367  n_components,
368  select,
369  integrate_to_one,
370  integral_Cinfty[dim - 1])
371  {}
372 
373 
374  template <int dim>
375  double
377  const unsigned int component) const
378  {
380  component == this->selected)
381  {
382  const double d = this->center.distance(p);
383  const double r = this->radius;
384  if (d >= r)
385  return 0.;
386  const double e = -r * r / (r * r - d * d);
387  return ((e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling);
388  }
389  return 0.;
390  }
391 
392 
393  template <int dim>
394  void
396  std::vector<double> & values,
397  const unsigned int component) const
398  {
399  Assert(values.size() == points.size(),
400  ExcDimensionMismatch(values.size(), points.size()));
401 
402  const double r = this->radius;
403 
405  component == this->selected)
406  for (unsigned int i = 0; i < values.size(); ++i)
407  {
408  const double d = this->center.distance(points[i]);
409  if (d >= r)
410  {
411  values[i] = 0.;
412  }
413  else
414  {
415  const double e = -r * r / (r * r - d * d);
416  values[i] =
417  (e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling;
418  }
419  }
420  else
421  std::fill(values.begin(), values.end(), 0.);
422  }
423 
424 
425  template <int dim>
426  void
428  const std::vector<Point<dim>> &points,
429  std::vector<Vector<double>> & values) const
430  {
431  Assert(values.size() == points.size(),
432  ExcDimensionMismatch(values.size(), points.size()));
433 
434  for (unsigned int k = 0; k < values.size(); ++k)
435  {
436  const double d = this->center.distance(points[k]);
437  const double r = this->radius;
438  double val = 0.;
439  if (d < this->radius)
440  {
441  const double e = -r * r / (r * r - d * d);
442  if (e > -50)
443  val = numbers::E * std::exp(e) * this->rescaling;
444  }
445 
447  values[k] = val;
448  else
449  {
450  values[k] = 0;
451  values[k](this->selected) = val;
452  }
453  }
454  }
455 
456 
457 
458  template <int dim>
461  const unsigned int) const
462  {
463  const double d = this->center.distance(p);
464  const double r = this->radius;
465  if (d >= r)
466  return Tensor<1, dim>();
467  const double e = -d * d / (r - d) / (r + d);
468  return ((e < -50) ? Point<dim>() :
469  (p - this->center) / d *
470  (-2.0 * r * r / std::pow(-r * r + d * d, 2.0) * d *
471  std::exp(e)) *
472  this->rescaling);
473  }
474 
475 
476 
477  template <int dim>
479  const Point<dim> p,
480  const unsigned int n_components,
481  const unsigned int select,
482  bool integrate_to_one)
483  : CutOffFunctionBase<dim>(r,
484  p,
485  n_components,
486  select,
487  integrate_to_one,
488  integral_C1[dim - 1])
489  {}
490 
491 
492  template <int dim>
493  double
495  const unsigned int component) const
496  {
498  component == this->selected)
499  {
500  const double d = this->center.distance(p);
501  const double r = this->radius;
502  if (d >= r)
503  return 0.;
504  return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
505  }
506  return 0.;
507  }
508 
509 
510  template <int dim>
511  void
512  CutOffFunctionC1<dim>::value_list(const std::vector<Point<dim>> &points,
513  std::vector<double> & values,
514  const unsigned int component) const
515  {
516  Assert(values.size() == points.size(),
517  ExcDimensionMismatch(values.size(), points.size()));
518 
519  const double r = this->radius;
520 
522  component == this->selected)
523  for (unsigned int i = 0; i < values.size(); ++i)
524  {
525  const double d = this->center.distance(points[i]);
526  if (d >= r)
527  {
528  values[i] = 0.;
529  }
530  else
531  {
532  values[i] =
533  .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
534  }
535  }
536  else
537  std::fill(values.begin(), values.end(), 0.);
538  }
539 
540 
541  template <int dim>
542  void
544  const std::vector<Point<dim>> &points,
545  std::vector<Vector<double>> & values) const
546  {
547  Assert(values.size() == points.size(),
548  ExcDimensionMismatch(values.size(), points.size()));
549 
550  for (unsigned int k = 0; k < values.size(); ++k)
551  {
552  const double d = this->center.distance(points[k]);
553  const double r = this->radius;
554  double val = 0.;
555  if (d < this->radius)
556  {
557  val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
558  }
559 
561  values[k] = val;
562  else
563  {
564  values[k] = 0;
565  values[k](this->selected) = val;
566  }
567  }
568  }
569 
570 
571 
572  template <int dim>
574  CutOffFunctionC1<dim>::gradient(const Point<dim> &p, const unsigned int) const
575  {
576  const double d = this->center.distance(p);
577  const double r = this->radius;
578  if (d >= r)
579  return Tensor<1, dim>();
580  return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) *
581  (p - this->center) / d * this->rescaling;
582  }
583 
584 
585  // explicit instantiations
586  template class CutOffFunctionBase<1>;
587  template class CutOffFunctionBase<2>;
588  template class CutOffFunctionBase<3>;
589 
590  template class CutOffFunctionLinfty<1>;
591  template class CutOffFunctionLinfty<2>;
592  template class CutOffFunctionLinfty<3>;
593 
594  template class CutOffFunctionW1<1>;
595  template class CutOffFunctionW1<2>;
596  template class CutOffFunctionW1<3>;
597 
598  template class CutOffFunctionCinfty<1>;
599  template class CutOffFunctionCinfty<2>;
600  template class CutOffFunctionCinfty<3>;
601 
602  template class CutOffFunctionC1<1>;
603  template class CutOffFunctionC1<2>;
604  template class CutOffFunctionC1<3>;
605 
606  template class CutOffFunctionTensorProduct<1>;
607  template class CutOffFunctionTensorProduct<2>;
608  template class CutOffFunctionTensorProduct<3>;
609 } // namespace Functions
610 
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:1691
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:1466
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:372
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: function_lib.h:989
Definition: cuda.h:32
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:371
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
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)