Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
function_lib_cutoff.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2001 - 2023 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
27namespace 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 {
124 Assert(initialized, ExcNotInitialized());
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 {
136 Assert(initialized, ExcNotInitialized());
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 {
149 Assert(initialized, ExcNotInitialized());
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 {
163 Assert(initialized, ExcNotInitialized());
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
177 //--------------------------------------------------------------------
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 {
227 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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
245 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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.;
270 if (this->selected == CutOffFunctionBase<dim>::no_component)
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 {
300 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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
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
321 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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.;
352 if (this->selected == CutOffFunctionBase<dim>::no_component)
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 {
384 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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
409 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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
451 if (this->selected == CutOffFunctionBase<dim>::no_component)
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) ?
474 Point<dim>() :
475 (p - this->center) / d *
476 (-2.0 * r * r / Utilities::fixed_power<2>(-r * r + d * d) * d *
477 std::exp(e)) *
478 this->rescaling);
479 }
480
481
482
483 template <int dim>
485 const Point<dim> p,
486 const unsigned int n_components,
487 const unsigned int select,
488 bool integrate_to_one)
489 : CutOffFunctionBase<dim>(r,
490 p,
491 n_components,
492 select,
493 integrate_to_one,
494 integral_C1[dim - 1])
495 {}
496
497
498 template <int dim>
499 double
501 const unsigned int component) const
502 {
503 if (this->selected == CutOffFunctionBase<dim>::no_component ||
504 component == this->selected)
505 {
506 const double d = this->center.distance(p);
507 const double r = this->radius;
508 if (d >= r)
509 return 0.;
510 return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
511 }
512 return 0.;
513 }
514
515
516 template <int dim>
517 void
519 std::vector<double> & values,
520 const unsigned int component) const
521 {
522 Assert(values.size() == points.size(),
523 ExcDimensionMismatch(values.size(), points.size()));
524
525 const double r = this->radius;
526
527 if (this->selected == CutOffFunctionBase<dim>::no_component ||
528 component == this->selected)
529 for (unsigned int i = 0; i < values.size(); ++i)
530 {
531 const double d = this->center.distance(points[i]);
532 if (d >= r)
533 {
534 values[i] = 0.;
535 }
536 else
537 {
538 values[i] =
539 .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
540 }
541 }
542 else
543 std::fill(values.begin(), values.end(), 0.);
544 }
545
546
547 template <int dim>
548 void
550 const std::vector<Point<dim>> &points,
551 std::vector<Vector<double>> & values) const
552 {
553 Assert(values.size() == points.size(),
554 ExcDimensionMismatch(values.size(), points.size()));
555
556 for (unsigned int k = 0; k < values.size(); ++k)
557 {
558 const double d = this->center.distance(points[k]);
559 const double r = this->radius;
560 double val = 0.;
561 if (d < this->radius)
562 {
563 val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
564 }
565
566 if (this->selected == CutOffFunctionBase<dim>::no_component)
567 values[k] = val;
568 else
569 {
570 values[k] = 0;
571 values[k](this->selected) = val;
572 }
573 }
574 }
575
576
577
578 template <int dim>
580 CutOffFunctionC1<dim>::gradient(const Point<dim> &p, const unsigned int) const
581 {
582 const double d = this->center.distance(p);
583 const double r = this->radius;
584 if (d >= r)
585 return Tensor<1, dim>();
586 return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) *
587 (p - this->center) / d * this->rescaling;
588 }
589
590
591 // explicit instantiations
592 template class CutOffFunctionBase<1>;
593 template class CutOffFunctionBase<2>;
594 template class CutOffFunctionBase<3>;
595
596 template class CutOffFunctionLinfty<1>;
597 template class CutOffFunctionLinfty<2>;
598 template class CutOffFunctionLinfty<3>;
599
600 template class CutOffFunctionW1<1>;
601 template class CutOffFunctionW1<2>;
602 template class CutOffFunctionW1<3>;
603
604 template class CutOffFunctionCinfty<1>;
605 template class CutOffFunctionCinfty<2>;
606 template class CutOffFunctionCinfty<3>;
607
608 template class CutOffFunctionC1<1>;
609 template class CutOffFunctionC1<2>;
610 template class CutOffFunctionC1<3>;
611
612 template class CutOffFunctionTensorProduct<1>;
613 template class CutOffFunctionTensorProduct<2>;
614 template class CutOffFunctionTensorProduct<3>;
615} // namespace Functions
616
virtual void set_radius(const double r)
virtual void set_center(const Point< dim > &p)
const Point< dim > & get_center() const
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)
virtual double value(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
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 Tensor< 1, dim > gradient(const Point< dim > &p, 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 void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, 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 Tensor< 1, dim > gradient(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
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 void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
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 void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
CutOffFunctionTensorProduct(double radius=1.0, 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)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void set_radius(const double radius) override
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
virtual void set_center(const Point< dim > &center) 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
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)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
Definition point.h:112
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > center
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
static constexpr double E
Definition numbers.h:234
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)