Reference documentation for deal.II version 9.3.3
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
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
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) ? 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 {
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
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 3 > center
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T fixed_power(const T t)
Definition: utilities.h:1081
static constexpr double E
Definition: numbers.h:206
static constexpr double PI
Definition: numbers.h:231
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)