Reference documentation for deal.II version Git 418c2197d0 2021-10-27 08:19:48 +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 
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 d = 0; d < dim; ++d)
166  {
167  ret[d] = base[d]->gradient(Point<1>(p[d]), component)[0];
168  for (unsigned int i = 0; i < dim; ++i)
169  if (i != d)
170  ret[d] *= base[i]->value(Point<1>(p[i]), component);
171  }
172  return ret;
173  }
174 
175 
176 
178  namespace
179  {
180  // Integral of CutOffFunctionLinfty in dimension 1, 2, and 3 when the radius
181  // is one
182  const double integral_Linfty[] = {2.0,
183  3.14159265358979323846264338328,
184  4.18879020478639098461685784437};
185 
186  // Integral of CutOffFunctionW1 in dimension 1, 2, and 3 when the radius
187  // is one
188  const double integral_W1[] = {1.0,
189  1.04719755119659774615421446109,
190  1.04719755119659774615421446109};
191 
192  // Integral of CutOffFunctionCinfty in dimension 1, 2, and 3 when the radius
193  // is one
194  const double integral_Cinfty[] = {1.20690032243787617533623799633,
195  1.26811216112759608094632335664,
196  1.1990039070192139033798473858};
197 
198  // Integral of CutOffFunctionC1 in dimension 1, 2, and 3 when the radius
199  // is one
200  const double integral_C1[] = {1.0,
201  0.93417655442731527615578663815,
202  0.821155557658032806157358815206};
203  } // namespace
204 
205 
206  template <int dim>
208  const double r,
209  const Point<dim> p,
210  const unsigned int n_components,
211  const unsigned int select,
212  const bool integrate_to_one)
213  : CutOffFunctionBase<dim>(r,
214  p,
215  n_components,
216  select,
217  integrate_to_one,
218  integral_Linfty[dim - 1])
219  {}
220 
221 
222  template <int dim>
223  double
225  const unsigned int component) const
226  {
228  component == this->selected)
229  return ((this->center.distance(p) < this->radius) ? this->rescaling : 0.);
230  return 0.;
231  }
232 
233 
234  template <int dim>
235  void
237  std::vector<double> & values,
238  const unsigned int component) const
239  {
240  Assert(values.size() == points.size(),
241  ExcDimensionMismatch(values.size(), points.size()));
242  AssertIndexRange(component, this->n_components);
243 
244 
246  component == this->selected)
247  for (unsigned int k = 0; k < values.size(); ++k)
248  values[k] = (this->center.distance(points[k]) < this->radius) ?
249  this->rescaling :
250  0.;
251  else
252  std::fill(values.begin(), values.end(), 0.);
253  }
254 
255 
256  template <int dim>
257  void
259  const std::vector<Point<dim>> &points,
260  std::vector<Vector<double>> & values) const
261  {
262  Assert(values.size() == points.size(),
263  ExcDimensionMismatch(values.size(), points.size()));
264 
265  for (unsigned int k = 0; k < values.size(); ++k)
266  {
267  const double val = (this->center.distance(points[k]) < this->radius) ?
268  this->rescaling :
269  0.;
271  values[k] = val;
272  else
273  {
274  values[k] = 0;
275  values[k](this->selected) = val;
276  }
277  }
278  }
279 
280  template <int dim>
282  const Point<dim> p,
283  const unsigned int n_components,
284  const unsigned int select,
285  const bool integrate_to_one)
286  : CutOffFunctionBase<dim>(r,
287  p,
288  n_components,
289  select,
290  integrate_to_one,
291  integral_W1[dim - 1])
292  {}
293 
294 
295  template <int dim>
296  double
298  const unsigned int component) const
299  {
301  component == this->selected)
302  {
303  const double d = this->center.distance(p);
304  return ((d < this->radius) ?
305  (this->radius - d) / this->radius * this->rescaling :
306  0.);
307  }
308  return 0.;
309  }
310 
311 
312  template <int dim>
313  void
314  CutOffFunctionW1<dim>::value_list(const std::vector<Point<dim>> &points,
315  std::vector<double> & values,
316  const unsigned int component) const
317  {
318  Assert(values.size() == points.size(),
319  ExcDimensionMismatch(values.size(), points.size()));
320 
322  component == this->selected)
323  for (unsigned int i = 0; i < values.size(); ++i)
324  {
325  const double d = this->center.distance(points[i]);
326  values[i] = ((d < this->radius) ?
327  (this->radius - d) / this->radius * this->rescaling :
328  0.);
329  }
330  else
331  std::fill(values.begin(), values.end(), 0.);
332  }
333 
334 
335 
336  template <int dim>
337  void
339  const std::vector<Point<dim>> &points,
340  std::vector<Vector<double>> & values) const
341  {
342  Assert(values.size() == points.size(),
343  ExcDimensionMismatch(values.size(), points.size()));
344 
345  for (unsigned int k = 0; k < values.size(); ++k)
346  {
347  const double d = this->center.distance(points[k]);
348  const double val =
349  (d < this->radius) ?
350  (this->radius - d) / this->radius * this->rescaling :
351  0.;
353  values[k] = val;
354  else
355  {
356  values[k] = 0;
357  values[k](this->selected) = val;
358  }
359  }
360  }
361 
362 
363  template <int dim>
365  const double r,
366  const Point<dim> p,
367  const unsigned int n_components,
368  const unsigned int select,
369  bool integrate_to_one)
370  : CutOffFunctionBase<dim>(r,
371  p,
372  n_components,
373  select,
374  integrate_to_one,
375  integral_Cinfty[dim - 1])
376  {}
377 
378 
379  template <int dim>
380  double
382  const unsigned int component) const
383  {
385  component == this->selected)
386  {
387  const double d = this->center.distance(p);
388  const double r = this->radius;
389  if (d >= r)
390  return 0.;
391  const double e = -r * r / (r * r - d * d);
392  return ((e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling);
393  }
394  return 0.;
395  }
396 
397 
398  template <int dim>
399  void
401  std::vector<double> & values,
402  const unsigned int component) const
403  {
404  Assert(values.size() == points.size(),
405  ExcDimensionMismatch(values.size(), points.size()));
406 
407  const double r = this->radius;
408 
410  component == this->selected)
411  for (unsigned int i = 0; i < values.size(); ++i)
412  {
413  const double d = this->center.distance(points[i]);
414  if (d >= r)
415  {
416  values[i] = 0.;
417  }
418  else
419  {
420  const double e = -r * r / (r * r - d * d);
421  values[i] =
422  (e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling;
423  }
424  }
425  else
426  std::fill(values.begin(), values.end(), 0.);
427  }
428 
429 
430  template <int dim>
431  void
433  const std::vector<Point<dim>> &points,
434  std::vector<Vector<double>> & values) const
435  {
436  Assert(values.size() == points.size(),
437  ExcDimensionMismatch(values.size(), points.size()));
438 
439  for (unsigned int k = 0; k < values.size(); ++k)
440  {
441  const double d = this->center.distance(points[k]);
442  const double r = this->radius;
443  double val = 0.;
444  if (d < this->radius)
445  {
446  const double e = -r * r / (r * r - d * d);
447  if (e > -50)
448  val = numbers::E * std::exp(e) * this->rescaling;
449  }
450 
452  values[k] = val;
453  else
454  {
455  values[k] = 0;
456  values[k](this->selected) = val;
457  }
458  }
459  }
460 
461 
462 
463  template <int dim>
466  const unsigned int) const
467  {
468  const double d = this->center.distance(p);
469  const double r = this->radius;
470  if (d >= r)
471  return Tensor<1, dim>();
472  const double e = -d * d / (r - d) / (r + d);
473  return ((e < -50) ? Point<dim>() :
474  (p - this->center) / d *
475  (-2.0 * r * r / std::pow(-r * r + d * d, 2.0) * d *
476  std::exp(e)) *
477  this->rescaling);
478  }
479 
480 
481 
482  template <int dim>
484  const Point<dim> p,
485  const unsigned int n_components,
486  const unsigned int select,
487  bool integrate_to_one)
488  : CutOffFunctionBase<dim>(r,
489  p,
490  n_components,
491  select,
492  integrate_to_one,
493  integral_C1[dim - 1])
494  {}
495 
496 
497  template <int dim>
498  double
500  const unsigned int component) const
501  {
503  component == this->selected)
504  {
505  const double d = this->center.distance(p);
506  const double r = this->radius;
507  if (d >= r)
508  return 0.;
509  return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
510  }
511  return 0.;
512  }
513 
514 
515  template <int dim>
516  void
517  CutOffFunctionC1<dim>::value_list(const std::vector<Point<dim>> &points,
518  std::vector<double> & values,
519  const unsigned int component) const
520  {
521  Assert(values.size() == points.size(),
522  ExcDimensionMismatch(values.size(), points.size()));
523 
524  const double r = this->radius;
525 
527  component == this->selected)
528  for (unsigned int i = 0; i < values.size(); ++i)
529  {
530  const double d = this->center.distance(points[i]);
531  if (d >= r)
532  {
533  values[i] = 0.;
534  }
535  else
536  {
537  values[i] =
538  .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
539  }
540  }
541  else
542  std::fill(values.begin(), values.end(), 0.);
543  }
544 
545 
546  template <int dim>
547  void
549  const std::vector<Point<dim>> &points,
550  std::vector<Vector<double>> & values) const
551  {
552  Assert(values.size() == points.size(),
553  ExcDimensionMismatch(values.size(), points.size()));
554 
555  for (unsigned int k = 0; k < values.size(); ++k)
556  {
557  const double d = this->center.distance(points[k]);
558  const double r = this->radius;
559  double val = 0.;
560  if (d < this->radius)
561  {
562  val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
563  }
564 
566  values[k] = val;
567  else
568  {
569  values[k] = 0;
570  values[k](this->selected) = val;
571  }
572  }
573  }
574 
575 
576 
577  template <int dim>
579  CutOffFunctionC1<dim>::gradient(const Point<dim> &p, const unsigned int) const
580  {
581  const double d = this->center.distance(p);
582  const double r = this->radius;
583  if (d >= r)
584  return Tensor<1, dim>();
585  return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) *
586  (p - this->center) / d * this->rescaling;
587  }
588 
589 
590  // explicit instantiations
591  template class CutOffFunctionBase<1>;
592  template class CutOffFunctionBase<2>;
593  template class CutOffFunctionBase<3>;
594 
595  template class CutOffFunctionLinfty<1>;
596  template class CutOffFunctionLinfty<2>;
597  template class CutOffFunctionLinfty<3>;
598 
599  template class CutOffFunctionW1<1>;
600  template class CutOffFunctionW1<2>;
601  template class CutOffFunctionW1<3>;
602 
603  template class CutOffFunctionCinfty<1>;
604  template class CutOffFunctionCinfty<2>;
605  template class CutOffFunctionCinfty<3>;
606 
607  template class CutOffFunctionC1<1>;
608  template class CutOffFunctionC1<2>;
609  template class CutOffFunctionC1<3>;
610 
611  template class CutOffFunctionTensorProduct<1>;
612  template class CutOffFunctionTensorProduct<2>;
613  template class CutOffFunctionTensorProduct<3>;
614 } // namespace Functions
615 
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:1720
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:1082
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:1461
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:401
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
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:400
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)