Reference documentation for deal.II version 9.6.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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16#include <deal.II/base/point.h>
17#include <deal.II/base/tensor.h>
18
19#include <deal.II/lac/vector.h>
20
21#include <cmath>
22
24
25
26namespace Functions
27{
28 template <int dim>
30 const double r,
31 const Point<dim> p,
32 const unsigned int n_components,
33 const unsigned int select,
34 const bool integrate_to_one,
35 const double unitary_integral_value)
36 : Function<dim>(n_components)
37 , center(p)
38 , radius(r)
39 , selected(select)
40 , integrate_to_one(integrate_to_one)
41 , unitary_integral_value(unitary_integral_value)
42 , rescaling(integrate_to_one ? 1. / (unitary_integral_value *
43 Utilities::fixed_power<dim>(radius)) :
44 1.0)
45 {
46 Assert(r > 0, ExcMessage("You must specify a radius > 0."));
47 }
48
49
50
51 template <int dim>
52 void
57
58
59
60 template <int dim>
61 const Point<dim> &
63 {
64 return center;
65 }
66
67
68
69 template <int dim>
70 void
72 {
73 radius = r;
74 Assert(r > 0, ExcMessage("You must specify a radius > 0."));
75 if (integrate_to_one)
76 rescaling =
77 1. / (unitary_integral_value * Utilities::fixed_power<dim>(radius));
78 else
79 rescaling = 1.0;
80 }
81
82
83
84 template <int dim>
85 double
87 {
88 return radius;
89 }
90
91
92
93 template <int dim>
94 bool
96 {
97 return integrate_to_one;
98 }
99
100
101
102 template <int dim>
104 double radius,
105 const Point<dim> &center,
106 const unsigned int n_components,
107 const unsigned int select,
108 const bool integrate_to_one)
109 : CutOffFunctionBase<dim>(radius,
110 center,
111 n_components,
112 select,
113 integrate_to_one)
114 , initialized(false)
115 {}
116
117
118
119 template <int dim>
120 void
122 {
123 Assert(initialized, ExcNotInitialized());
124 for (unsigned int i = 0; i < dim; ++i)
125 base[i]->set_center(Point<1>(p[i]));
127 }
128
129
130
131 template <int dim>
132 void
134 {
135 Assert(initialized, ExcNotInitialized());
136 for (unsigned int i = 0; i < dim; ++i)
137 base[i]->set_radius(r);
139 }
140
141
142
143 template <int dim>
144 double
146 const unsigned int component) const
147 {
148 Assert(initialized, ExcNotInitialized());
149 double ret = 1.0;
150 for (unsigned int i = 0; i < dim; ++i)
151 ret *= base[i]->value(Point<1>(p[i]), component);
152 return ret;
153 }
154
155
156
157 template <int dim>
160 const unsigned int component) const
161 {
162 Assert(initialized, ExcNotInitialized());
163 Tensor<1, dim> ret;
164 for (unsigned int d = 0; d < dim; ++d)
165 {
166 ret[d] = base[d]->gradient(Point<1>(p[d]), component)[0];
167 for (unsigned int i = 0; i < dim; ++i)
168 if (i != d)
169 ret[d] *= base[i]->value(Point<1>(p[i]), component);
170 }
171 return ret;
172 }
173
174
175
176 //--------------------------------------------------------------------
177 namespace
178 {
179 // Integral of CutOffFunctionLinfty in dimension 1, 2, and 3 when the radius
180 // is one
181 const double integral_Linfty[] = {2.0,
182 3.14159265358979323846264338328,
183 4.18879020478639098461685784437};
184
185 // Integral of CutOffFunctionW1 in dimension 1, 2, and 3 when the radius
186 // is one
187 const double integral_W1[] = {1.0,
188 1.04719755119659774615421446109,
189 1.04719755119659774615421446109};
190
191 // Integral of CutOffFunctionCinfty in dimension 1, 2, and 3 when the radius
192 // is one
193 const double integral_Cinfty[] = {1.20690032243787617533623799633,
194 1.26811216112759608094632335664,
195 1.1990039070192139033798473858};
196
197 // Integral of CutOffFunctionC1 in dimension 1, 2, and 3 when the radius
198 // is one
199 const double integral_C1[] = {1.0,
200 0.93417655442731527615578663815,
201 0.821155557658032806157358815206};
202 } // namespace
203
204
205 template <int dim>
207 const double r,
208 const Point<dim> p,
209 const unsigned int n_components,
210 const unsigned int select,
211 const bool integrate_to_one)
212 : CutOffFunctionBase<dim>(r,
213 p,
214 n_components,
215 select,
216 integrate_to_one,
217 integral_Linfty[dim - 1])
218 {}
219
220
221 template <int dim>
222 double
224 const unsigned int component) const
225 {
226 if (this->selected == CutOffFunctionBase<dim>::no_component ||
227 component == this->selected)
228 return ((this->center.distance(p) < this->radius) ? this->rescaling : 0.);
229 return 0.;
230 }
231
232
233 template <int dim>
234 void
236 std::vector<double> &values,
237 const unsigned int component) const
238 {
239 Assert(values.size() == points.size(),
240 ExcDimensionMismatch(values.size(), points.size()));
241 AssertIndexRange(component, this->n_components);
242
243
244 if (this->selected == CutOffFunctionBase<dim>::no_component ||
245 component == this->selected)
246 for (unsigned int k = 0; k < values.size(); ++k)
247 values[k] = (this->center.distance(points[k]) < this->radius) ?
248 this->rescaling :
249 0.;
250 else
251 std::fill(values.begin(), values.end(), 0.);
252 }
253
254
255 template <int dim>
256 void
258 const std::vector<Point<dim>> &points,
259 std::vector<Vector<double>> &values) const
260 {
261 Assert(values.size() == points.size(),
262 ExcDimensionMismatch(values.size(), points.size()));
263
264 for (unsigned int k = 0; k < values.size(); ++k)
265 {
266 const double val = (this->center.distance(points[k]) < this->radius) ?
267 this->rescaling :
268 0.;
269 if (this->selected == CutOffFunctionBase<dim>::no_component)
270 values[k] = val;
271 else
272 {
273 values[k] = 0;
274 values[k](this->selected) = val;
275 }
276 }
277 }
278
279 template <int dim>
281 const Point<dim> p,
282 const unsigned int n_components,
283 const unsigned int select,
284 const bool integrate_to_one)
285 : CutOffFunctionBase<dim>(r,
286 p,
287 n_components,
288 select,
289 integrate_to_one,
290 integral_W1[dim - 1])
291 {}
292
293
294 template <int dim>
295 double
297 const unsigned int component) const
298 {
299 if (this->selected == CutOffFunctionBase<dim>::no_component ||
300 component == this->selected)
301 {
302 const double d = this->center.distance(p);
303 return ((d < this->radius) ?
304 (this->radius - d) / this->radius * this->rescaling :
305 0.);
306 }
307 return 0.;
308 }
309
310
311 template <int dim>
312 void
314 std::vector<double> &values,
315 const unsigned int component) const
316 {
317 Assert(values.size() == points.size(),
318 ExcDimensionMismatch(values.size(), points.size()));
319
320 if (this->selected == CutOffFunctionBase<dim>::no_component ||
321 component == this->selected)
322 for (unsigned int i = 0; i < values.size(); ++i)
323 {
324 const double d = this->center.distance(points[i]);
325 values[i] = ((d < this->radius) ?
326 (this->radius - d) / this->radius * this->rescaling :
327 0.);
328 }
329 else
330 std::fill(values.begin(), values.end(), 0.);
331 }
332
333
334
335 template <int dim>
336 void
338 const std::vector<Point<dim>> &points,
339 std::vector<Vector<double>> &values) const
340 {
341 Assert(values.size() == points.size(),
342 ExcDimensionMismatch(values.size(), points.size()));
343
344 for (unsigned int k = 0; k < values.size(); ++k)
345 {
346 const double d = this->center.distance(points[k]);
347 const double val =
348 (d < this->radius) ?
349 (this->radius - d) / this->radius * this->rescaling :
350 0.;
351 if (this->selected == CutOffFunctionBase<dim>::no_component)
352 values[k] = val;
353 else
354 {
355 values[k] = 0;
356 values[k](this->selected) = val;
357 }
358 }
359 }
360
361
362 template <int dim>
364 const double r,
365 const Point<dim> p,
366 const unsigned int n_components,
367 const unsigned int select,
368 bool integrate_to_one)
369 : CutOffFunctionBase<dim>(r,
370 p,
371 n_components,
372 select,
373 integrate_to_one,
374 integral_Cinfty[dim - 1])
375 {}
376
377
378 template <int dim>
379 double
381 const unsigned int component) const
382 {
383 if (this->selected == CutOffFunctionBase<dim>::no_component ||
384 component == this->selected)
385 {
386 const double d = this->center.distance(p);
387 const double r = this->radius;
388 if (d >= r)
389 return 0.;
390 const double e = -r * r / (r * r - d * d);
391 return ((e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling);
392 }
393 return 0.;
394 }
395
396
397 template <int dim>
398 void
400 std::vector<double> &values,
401 const unsigned int component) const
402 {
403 Assert(values.size() == points.size(),
404 ExcDimensionMismatch(values.size(), points.size()));
405
406 const double r = this->radius;
407
408 if (this->selected == CutOffFunctionBase<dim>::no_component ||
409 component == this->selected)
410 for (unsigned int i = 0; i < values.size(); ++i)
411 {
412 const double d = this->center.distance(points[i]);
413 if (d >= r)
414 {
415 values[i] = 0.;
416 }
417 else
418 {
419 const double e = -r * r / (r * r - d * d);
420 values[i] =
421 (e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling;
422 }
423 }
424 else
425 std::fill(values.begin(), values.end(), 0.);
426 }
427
428
429 template <int dim>
430 void
432 const std::vector<Point<dim>> &points,
433 std::vector<Vector<double>> &values) const
434 {
435 Assert(values.size() == points.size(),
436 ExcDimensionMismatch(values.size(), points.size()));
437
438 for (unsigned int k = 0; k < values.size(); ++k)
439 {
440 const double d = this->center.distance(points[k]);
441 const double r = this->radius;
442 double val = 0.;
443 if (d < this->radius)
444 {
445 const double e = -r * r / (r * r - d * d);
446 if (e > -50)
447 val = numbers::E * std::exp(e) * this->rescaling;
448 }
449
450 if (this->selected == CutOffFunctionBase<dim>::no_component)
451 values[k] = val;
452 else
453 {
454 values[k] = 0;
455 values[k](this->selected) = val;
456 }
457 }
458 }
459
460
461
462 template <int dim>
465 const unsigned int) const
466 {
467 const double d = this->center.distance(p);
468 const double r = this->radius;
469 if (d >= r)
470 return Tensor<1, dim>();
471 const double e = -d * d / (r - d) / (r + d);
472 return ((e < -50) ?
473 Point<dim>() :
474 (p - this->center) / d *
475 (-2.0 * r * r / Utilities::fixed_power<2>(-r * r + d * d) * 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 {
502 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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
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
526 if (this->selected == CutOffFunctionBase<dim>::no_component ||
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
565 if (this->selected == CutOffFunctionBase<dim>::no_component)
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 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:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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)
constexpr T fixed_power(const T t)
Definition utilities.h:942
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 > &)