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
flow_function.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2007 - 2024 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
23
25
26
27namespace Functions
28{
29 template <int dim>
31 : Function<dim>(dim + 1)
32 , mean_pressure(0)
33 , aux_values(dim + 1)
34 , aux_gradients(dim + 1)
35 {}
36
37
38
39 template <int dim>
40 void
42 {
43 mean_pressure = p;
44 }
45
46
47 template <int dim>
48 void
50 const std::vector<Point<dim>> &points,
51 std::vector<Vector<double>> &values) const
52 {
53 const unsigned int n_points = points.size();
54 Assert(values.size() == n_points,
55 ExcDimensionMismatch(values.size(), n_points));
56
57 // guard access to the aux_*
58 // variables in multithread mode
59 std::lock_guard<std::mutex> lock(mutex);
60
61 for (unsigned int d = 0; d < dim + 1; ++d)
62 aux_values[d].resize(n_points);
63 vector_values(points, aux_values);
64
65 for (unsigned int k = 0; k < n_points; ++k)
66 {
67 Assert(values[k].size() == dim + 1,
68 ExcDimensionMismatch(values[k].size(), dim + 1));
69 for (unsigned int d = 0; d < dim + 1; ++d)
70 values[k](d) = aux_values[d][k];
71 }
72 }
73
74
75 template <int dim>
76 void
78 Vector<double> &value) const
79 {
80 Assert(value.size() == dim + 1,
81 ExcDimensionMismatch(value.size(), dim + 1));
82
83 const unsigned int n_points = 1;
84 std::vector<Point<dim>> points(1);
85 points[0] = point;
86
87 // guard access to the aux_*
88 // variables in multithread mode
89 std::lock_guard<std::mutex> lock(mutex);
90
91 for (unsigned int d = 0; d < dim + 1; ++d)
92 aux_values[d].resize(n_points);
93 vector_values(points, aux_values);
94
95 for (unsigned int d = 0; d < dim + 1; ++d)
96 value(d) = aux_values[d][0];
97 }
98
99
100 template <int dim>
101 double
103 const unsigned int comp) const
104 {
105 AssertIndexRange(comp, dim + 1);
106 const unsigned int n_points = 1;
107 std::vector<Point<dim>> points(1);
108 points[0] = point;
109
110 // guard access to the aux_*
111 // variables in multithread mode
112 std::lock_guard<std::mutex> lock(mutex);
114 for (unsigned int d = 0; d < dim + 1; ++d)
115 aux_values[d].resize(n_points);
116 vector_values(points, aux_values);
117
118 return aux_values[comp][0];
119 }
121
122 template <int dim>
123 void
125 const std::vector<Point<dim>> &points,
126 std::vector<std::vector<Tensor<1, dim>>> &values) const
127 {
128 const unsigned int n_points = points.size();
129 Assert(values.size() == n_points,
130 ExcDimensionMismatch(values.size(), n_points));
131
132 // guard access to the aux_*
133 // variables in multithread mode
134 std::lock_guard<std::mutex> lock(mutex);
135
136 for (unsigned int d = 0; d < dim + 1; ++d)
137 aux_gradients[d].resize(n_points);
138 vector_gradients(points, aux_gradients);
139
140 for (unsigned int k = 0; k < n_points; ++k)
141 {
142 Assert(values[k].size() == dim + 1,
143 ExcDimensionMismatch(values[k].size(), dim + 1));
144 for (unsigned int d = 0; d < dim + 1; ++d)
145 values[k][d] = aux_gradients[d][k];
146 }
147 }
148
149
150 template <int dim>
151 void
153 const std::vector<Point<dim>> &points,
154 std::vector<Vector<double>> &values) const
155 {
156 const unsigned int n_points = points.size();
157 Assert(values.size() == n_points,
158 ExcDimensionMismatch(values.size(), n_points));
159
160 // guard access to the aux_*
161 // variables in multithread mode
162 std::lock_guard<std::mutex> lock(mutex);
163
164 for (unsigned int d = 0; d < dim + 1; ++d)
165 aux_values[d].resize(n_points);
166 vector_laplacians(points, aux_values);
167
168 for (unsigned int k = 0; k < n_points; ++k)
169 {
170 Assert(values[k].size() == dim + 1,
171 ExcDimensionMismatch(values[k].size(), dim + 1));
172 for (unsigned int d = 0; d < dim + 1; ++d)
173 values[k](d) = aux_values[d][k];
174 }
175 }
176
177
178 template <int dim>
179 std::size_t
181 {
183 return 0;
184 }
185
186
187 //----------------------------------------------------------------------//
188
189 template <int dim>
190 PoisseuilleFlow<dim>::PoisseuilleFlow(const double r, const double Re)
191 : inv_sqr_radius(1 / r / r)
192 , Reynolds(Re)
193 {
194 Assert(Reynolds != 0., ExcMessage("Reynolds number cannot be zero"));
195 }
196
197
198
199 template <int dim>
200 void
202 const std::vector<Point<dim>> &points,
203 std::vector<std::vector<double>> &values) const
204 {
205 const unsigned int n = points.size();
206
207 Assert(values.size() == dim + 1,
208 ExcDimensionMismatch(values.size(), dim + 1));
209 for (unsigned int d = 0; d < dim + 1; ++d)
210 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
211
212 for (unsigned int k = 0; k < n; ++k)
213 {
214 const Point<dim> &p = points[k];
215 // First, compute the square of the distance to the x-axis divided by
216 // the radius.
217 double r2 = 0;
218 for (unsigned int d = 1; d < dim; ++d)
219 r2 += p[d] * p[d];
220 r2 *= inv_sqr_radius;
221
222 // x-velocity
223 values[0][k] = 1. - r2;
224 // other velocities
225 for (unsigned int d = 1; d < dim; ++d)
226 values[d][k] = 0.;
227 // pressure
228 values[dim][k] = -2 * (dim - 1) * inv_sqr_radius * p[0] / Reynolds +
229 this->mean_pressure;
230 }
231 }
232
233
234
235 template <int dim>
236 void
238 const std::vector<Point<dim>> &points,
239 std::vector<std::vector<Tensor<1, dim>>> &values) const
240 {
241 const unsigned int n = points.size();
242
243 Assert(values.size() == dim + 1,
244 ExcDimensionMismatch(values.size(), dim + 1));
245 for (unsigned int d = 0; d < dim + 1; ++d)
246 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
247
248 for (unsigned int k = 0; k < n; ++k)
249 {
250 const Point<dim> &p = points[k];
251 // x-velocity
252 values[0][k][0] = 0.;
253 for (unsigned int d = 1; d < dim; ++d)
254 values[0][k][d] = -2. * p[d] * inv_sqr_radius;
255 // other velocities
256 for (unsigned int d = 1; d < dim; ++d)
257 values[d][k] = 0.;
258 // pressure
259 values[dim][k][0] = -2 * (dim - 1) * inv_sqr_radius / Reynolds;
260 for (unsigned int d = 1; d < dim; ++d)
261 values[dim][k][d] = 0.;
262 }
263 }
264
265
266
267 template <int dim>
268 void
270 const std::vector<Point<dim>> &points,
271 std::vector<std::vector<double>> &values) const
272 {
273 const unsigned int n = points.size();
274 (void)n;
275 Assert(values.size() == dim + 1,
276 ExcDimensionMismatch(values.size(), dim + 1));
277 for (unsigned int d = 0; d < dim + 1; ++d)
278 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
279
280 for (auto &point_values : values)
281 std::fill(point_values.begin(), point_values.end(), 0.);
282 }
283
284 //----------------------------------------------------------------------//
285
286 template <int dim>
287 StokesCosine<dim>::StokesCosine(const double nu, const double r)
288 : viscosity(nu)
289 , reaction(r)
290 {}
291
292
293
294 template <int dim>
295 void
296 StokesCosine<dim>::set_parameters(const double nu, const double r)
297 {
298 viscosity = nu;
299 reaction = r;
300 }
301
302
303 template <int dim>
304 void
306 const std::vector<Point<dim>> &points,
307 std::vector<std::vector<double>> &values) const
308 {
309 unsigned int n = points.size();
310
311 Assert(values.size() == dim + 1,
312 ExcDimensionMismatch(values.size(), dim + 1));
313 for (unsigned int d = 0; d < dim + 1; ++d)
314 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
315
316 for (unsigned int k = 0; k < n; ++k)
317 {
318 const Point<dim> &p = points[k];
319 const double x = numbers::PI / 2. * p[0];
320 const double y = numbers::PI / 2. * p[1];
321 const double cx = std::cos(x);
322 const double cy = std::cos(y);
323 const double sx = std::sin(x);
324 const double sy = std::sin(y);
325
326 if (dim == 2)
327 {
328 values[0][k] = cx * cx * cy * sy;
329 values[1][k] = -cx * sx * cy * cy;
330 values[2][k] = cx * sx * cy * sy + this->mean_pressure;
331 }
332 else if (dim == 3)
333 {
334 const double z = numbers::PI / 2. * p[2];
335 const double cz = std::cos(z);
336 const double sz = std::sin(z);
337
338 values[0][k] = cx * cx * cy * sy * cz * sz;
339 values[1][k] = cx * sx * cy * cy * cz * sz;
340 values[2][k] = -2. * cx * sx * cy * sy * cz * cz;
341 values[3][k] = cx * sx * cy * sy * cz * sz + this->mean_pressure;
342 }
343 else
344 {
346 }
347 }
348 }
349
350
351
352 template <int dim>
353 void
355 const std::vector<Point<dim>> &points,
356 std::vector<std::vector<Tensor<1, dim>>> &values) const
357 {
358 unsigned int n = points.size();
359
360 Assert(values.size() == dim + 1,
361 ExcDimensionMismatch(values.size(), dim + 1));
362 for (unsigned int d = 0; d < dim + 1; ++d)
363 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
364
365 for (unsigned int k = 0; k < n; ++k)
366 {
367 const Point<dim> &p = points[k];
368 const double x = numbers::PI / 2. * p[0];
369 const double y = numbers::PI / 2. * p[1];
370 const double c2x = std::cos(2 * x);
371 const double c2y = std::cos(2 * y);
372 const double s2x = std::sin(2 * x);
373 const double s2y = std::sin(2 * y);
374 const double cx2 = .5 + .5 * c2x; // cos^2 x
375 const double cy2 = .5 + .5 * c2y; // cos^2 y
376
377 if (dim == 2)
378 {
379 values[0][k][0] = -.25 * numbers::PI * s2x * s2y;
380 values[0][k][1] = .5 * numbers::PI * cx2 * c2y;
381 values[1][k][0] = -.5 * numbers::PI * c2x * cy2;
382 values[1][k][1] = .25 * numbers::PI * s2x * s2y;
383 values[2][k][0] = .25 * numbers::PI * c2x * s2y;
384 values[2][k][1] = .25 * numbers::PI * s2x * c2y;
385 }
386 else if (dim == 3)
387 {
388 const double z = numbers::PI / 2. * p[2];
389 const double c2z = std::cos(2 * z);
390 const double s2z = std::sin(2 * z);
391 const double cz2 = .5 + .5 * c2z; // cos^2 z
392
393 values[0][k][0] = -.125 * numbers::PI * s2x * s2y * s2z;
394 values[0][k][1] = .25 * numbers::PI * cx2 * c2y * s2z;
395 values[0][k][2] = .25 * numbers::PI * cx2 * s2y * c2z;
396
397 values[1][k][0] = .25 * numbers::PI * c2x * cy2 * s2z;
398 values[1][k][1] = -.125 * numbers::PI * s2x * s2y * s2z;
399 values[1][k][2] = .25 * numbers::PI * s2x * cy2 * c2z;
400
401 values[2][k][0] = -.5 * numbers::PI * c2x * s2y * cz2;
402 values[2][k][1] = -.5 * numbers::PI * s2x * c2y * cz2;
403 values[2][k][2] = .25 * numbers::PI * s2x * s2y * s2z;
404
405 values[3][k][0] = .125 * numbers::PI * c2x * s2y * s2z;
406 values[3][k][1] = .125 * numbers::PI * s2x * c2y * s2z;
407 values[3][k][2] = .125 * numbers::PI * s2x * s2y * c2z;
408 }
409 else
410 {
412 }
413 }
414 }
415
416
417
418 template <int dim>
419 void
421 const std::vector<Point<dim>> &points,
422 std::vector<std::vector<double>> &values) const
423 {
424 unsigned int n = points.size();
425
426 Assert(values.size() == dim + 1,
427 ExcDimensionMismatch(values.size(), dim + 1));
428 for (unsigned int d = 0; d < dim + 1; ++d)
429 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
430
431 if (reaction != 0.)
432 {
433 vector_values(points, values);
434 for (unsigned int d = 0; d < dim; ++d)
435 for (double &point_value : values[d])
436 point_value *= -reaction;
437 }
438 else
439 {
440 for (unsigned int d = 0; d < dim; ++d)
441 std::fill(values[d].begin(), values[d].end(), 0.);
442 }
443
444
445 for (unsigned int k = 0; k < n; ++k)
446 {
447 const Point<dim> &p = points[k];
448 const double x = numbers::PI / 2. * p[0];
449 const double y = numbers::PI / 2. * p[1];
450 const double c2x = std::cos(2 * x);
451 const double c2y = std::cos(2 * y);
452 const double s2x = std::sin(2 * x);
453 const double s2y = std::sin(2 * y);
454 const double pi2 = .25 * numbers::PI * numbers::PI;
455
456 if (dim == 2)
457 {
458 values[0][k] += -viscosity * pi2 * (1. + 2. * c2x) * s2y -
459 numbers::PI / 4. * c2x * s2y;
460 values[1][k] += viscosity * pi2 * s2x * (1. + 2. * c2y) -
461 numbers::PI / 4. * s2x * c2y;
462 values[2][k] = 0.;
463 }
464 else if (dim == 3)
465 {
466 const double z = numbers::PI * p[2];
467 const double c2z = std::cos(2 * z);
468 const double s2z = std::sin(2 * z);
469
470 values[0][k] +=
471 -.5 * viscosity * pi2 * (1. + 2. * c2x) * s2y * s2z -
472 numbers::PI / 8. * c2x * s2y * s2z;
473 values[1][k] += .5 * viscosity * pi2 * s2x * (1. + 2. * c2y) * s2z -
474 numbers::PI / 8. * s2x * c2y * s2z;
475 values[2][k] +=
476 -.5 * viscosity * pi2 * s2x * s2y * (1. + 2. * c2z) -
477 numbers::PI / 8. * s2x * s2y * c2z;
478 values[3][k] = 0.;
479 }
480 else
481 {
483 }
484 }
485 }
486
487
488 //----------------------------------------------------------------------//
489
490 const double StokesLSingularity::lambda = 0.54448373678246;
491
493 : omega(3. / 2. * numbers::PI)
494 , coslo(std::cos(lambda * omega))
495 , lp(1. + lambda)
496 , lm(1. - lambda)
497 {}
498
499
500 inline double
501 StokesLSingularity::Psi(double phi) const
502 {
503 return coslo * (std::sin(lp * phi) / lp - std::sin(lm * phi) / lm) -
504 std::cos(lp * phi) + std::cos(lm * phi);
505 }
506
507
508 inline double
510 {
511 return coslo * (std::cos(lp * phi) - std::cos(lm * phi)) +
512 lp * std::sin(lp * phi) - lm * std::sin(lm * phi);
513 }
514
515
516 inline double
518 {
519 return coslo * (lm * std::sin(lm * phi) - lp * std::sin(lp * phi)) +
520 lp * lp * std::cos(lp * phi) - lm * lm * std::cos(lm * phi);
521 }
522
523
524 inline double
526 {
527 return coslo *
528 (lm * lm * std::cos(lm * phi) - lp * lp * std::cos(lp * phi)) +
529 lm * lm * lm * std::sin(lm * phi) -
530 lp * lp * lp * std::sin(lp * phi);
531 }
532
533
534 inline double
536 {
537 return coslo * (lp * lp * lp * std::sin(lp * phi) -
538 lm * lm * lm * std::sin(lm * phi)) +
539 lm * lm * lm * lm * std::cos(lm * phi) -
540 lp * lp * lp * lp * std::cos(lp * phi);
541 }
542
543
544 void
546 const std::vector<Point<2>> &points,
547 std::vector<std::vector<double>> &values) const
548 {
549 unsigned int n = points.size();
550
551 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
552 for (unsigned int d = 0; d < 2 + 1; ++d)
553 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
554
555 for (unsigned int k = 0; k < n; ++k)
556 {
557 const Point<2> &p = points[k];
558 const double x = p[0];
559 const double y = p[1];
560
561 if ((x < 0) || (y < 0))
562 {
563 const double phi = std::atan2(y, -x) + numbers::PI;
564 const double r2 = x * x + y * y;
565 const double rl = std::pow(r2, lambda / 2.);
566 const double rl1 = std::pow(r2, lambda / 2. - .5);
567 values[0][k] =
568 rl * (lp * std::sin(phi) * Psi(phi) + std::cos(phi) * Psi_1(phi));
569 values[1][k] =
570 rl * (lp * std::cos(phi) * Psi(phi) - std::sin(phi) * Psi_1(phi));
571 values[2][k] = -rl1 * (lp * lp * Psi_1(phi) + Psi_3(phi)) / lm +
572 this->mean_pressure;
573 }
574 else
575 {
576 for (unsigned int d = 0; d < 3; ++d)
577 values[d][k] = 0.;
578 }
579 }
580 }
581
582
583
584 void
586 const std::vector<Point<2>> &points,
587 std::vector<std::vector<Tensor<1, 2>>> &values) const
588 {
589 unsigned int n = points.size();
590
591 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
592 for (unsigned int d = 0; d < 2 + 1; ++d)
593 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
594
595 for (unsigned int k = 0; k < n; ++k)
596 {
597 const Point<2> &p = points[k];
598 const double x = p[0];
599 const double y = p[1];
600
601 if ((x < 0) || (y < 0))
602 {
603 const double phi = std::atan2(y, -x) + numbers::PI;
604 const double r2 = x * x + y * y;
605 const double r = std::sqrt(r2);
606 const double rl = std::pow(r2, lambda / 2.);
607 const double rl1 = std::pow(r2, lambda / 2. - .5);
608 const double rl2 = std::pow(r2, lambda / 2. - 1.);
609 const double psi = Psi(phi);
610 const double psi1 = Psi_1(phi);
611 const double psi2 = Psi_2(phi);
612 const double cosp = std::cos(phi);
613 const double sinp = std::sin(phi);
614
615 // Derivatives of u with respect to r, phi
616 const double udr = lambda * rl1 * (lp * sinp * psi + cosp * psi1);
617 const double udp = rl * (lp * cosp * psi + lp * sinp * psi1 -
618 sinp * psi1 + cosp * psi2);
619 // Derivatives of v with respect to r, phi
620 const double vdr = lambda * rl1 * (lp * cosp * psi - sinp * psi1);
621 const double vdp = rl * (lp * (cosp * psi1 - sinp * psi) -
622 cosp * psi1 - sinp * psi2);
623 // Derivatives of p with respect to r, phi
624 const double pdr =
625 -(lambda - 1.) * rl2 * (lp * lp * psi1 + Psi_3(phi)) / lm;
626 const double pdp = -rl1 * (lp * lp * psi2 + Psi_4(phi)) / lm;
627 values[0][k][0] = cosp * udr - sinp / r * udp;
628 values[0][k][1] = -sinp * udr - cosp / r * udp;
629 values[1][k][0] = cosp * vdr - sinp / r * vdp;
630 values[1][k][1] = -sinp * vdr - cosp / r * vdp;
631 values[2][k][0] = cosp * pdr - sinp / r * pdp;
632 values[2][k][1] = -sinp * pdr - cosp / r * pdp;
633 }
634 else
635 {
636 for (unsigned int d = 0; d < 3; ++d)
637 values[d][k] = 0.;
638 }
639 }
640 }
641
642
643
644 void
646 const std::vector<Point<2>> &points,
647 std::vector<std::vector<double>> &values) const
648 {
649 unsigned int n = points.size();
650 (void)n;
651 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
652 for (unsigned int d = 0; d < 2 + 1; ++d)
653 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
654
655 for (auto &point_values : values)
656 std::fill(point_values.begin(), point_values.end(), 0.);
657 }
658
659
660 //----------------------------------------------------------------------//
661
662 Kovasznay::Kovasznay(double Re, bool stokes)
663 : Reynolds(Re)
664 , stokes(stokes)
665 {
666 long double r2 = Reynolds / 2.;
667 long double b = 4 * numbers::PI * numbers::PI;
668 long double l = -b / (r2 + std::sqrt(r2 * r2 + b));
669 lbda = l;
670 // mean pressure for a domain
671 // spreading from -.5 to 1.5 in
672 // x-direction
673 p_average = 1 / (8 * l) * (std::exp(3. * l) - std::exp(-l));
674 }
675
676
677
678 void
679 Kovasznay::vector_values(const std::vector<Point<2>> &points,
680 std::vector<std::vector<double>> &values) const
681 {
682 unsigned int n = points.size();
683
684 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
685 for (unsigned int d = 0; d < 2 + 1; ++d)
686 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
687
688 for (unsigned int k = 0; k < n; ++k)
689 {
690 const Point<2> &p = points[k];
691 const double x = p[0];
692 const double y = 2. * numbers::PI * p[1];
693 const double elx = std::exp(lbda * x);
694
695 values[0][k] = 1. - elx * std::cos(y);
696 values[1][k] = .5 / numbers::PI * lbda * elx * std::sin(y);
697 values[2][k] = -.5 * elx * elx + p_average + this->mean_pressure;
698 }
699 }
700
701
702 void
704 const std::vector<Point<2>> &points,
705 std::vector<std::vector<Tensor<1, 2>>> &gradients) const
706 {
707 unsigned int n = points.size();
708
709 Assert(gradients.size() == 3, ExcDimensionMismatch(gradients.size(), 3));
710 Assert(gradients[0].size() == n,
711 ExcDimensionMismatch(gradients[0].size(), n));
712
713 for (unsigned int i = 0; i < n; ++i)
714 {
715 const double x = points[i][0];
716 const double y = points[i][1];
717
718 const double elx = std::exp(lbda * x);
719 const double cy = std::cos(2 * numbers::PI * y);
720 const double sy = std::sin(2 * numbers::PI * y);
721
722 // u
723 gradients[0][i][0] = -lbda * elx * cy;
724 gradients[0][i][1] = 2. * numbers::PI * elx * sy;
725 gradients[1][i][0] = lbda * lbda / (2 * numbers::PI) * elx * sy;
726 gradients[1][i][1] = lbda * elx * cy;
727 // p
728 gradients[2][i][0] = -lbda * elx * elx;
729 gradients[2][i][1] = 0.;
730 }
731 }
732
733
734
735 void
736 Kovasznay::vector_laplacians(const std::vector<Point<2>> &points,
737 std::vector<std::vector<double>> &values) const
738 {
739 unsigned int n = points.size();
740 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
741 for (unsigned int d = 0; d < 2 + 1; ++d)
742 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
743
744 if (stokes)
745 {
746 const double zp = 2. * numbers::PI;
747 for (unsigned int k = 0; k < n; ++k)
748 {
749 const Point<2> &p = points[k];
750 const double x = p[0];
751 const double y = zp * p[1];
752 const double elx = std::exp(lbda * x);
753 const double u = 1. - elx * std::cos(y);
754 const double ux = -lbda * elx * std::cos(y);
755 const double uy = elx * zp * std::sin(y);
756 const double v = lbda / zp * elx * std::sin(y);
757 const double vx = lbda * lbda / zp * elx * std::sin(y);
758 const double vy = zp * lbda / zp * elx * std::cos(y);
759
760 values[0][k] = u * ux + v * uy;
761 values[1][k] = u * vx + v * vy;
762 values[2][k] = 0.;
763 }
764 }
765 else
766 {
767 for (auto &point_values : values)
768 std::fill(point_values.begin(), point_values.end(), 0.);
769 }
770 }
771
772 double
774 {
775 return lbda;
776 }
777
778
779
780 template class FlowFunction<2>;
781 template class FlowFunction<3>;
782 template class PoisseuilleFlow<2>;
783 template class PoisseuilleFlow<3>;
784 template class StokesCosine<2>;
785 template class StokesCosine<3>;
786} // namespace Functions
787
788
789
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
void pressure_adjustment(double p)
virtual void vector_gradient_list(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const override
virtual void vector_laplacian_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const override
virtual std::size_t memory_consumption() const override
virtual void vector_value(const Point< dim > &points, Vector< double > &value) const override
virtual double value(const Point< dim > &points, const unsigned int component) const override
Kovasznay(const double Re, bool Stokes=false)
virtual void vector_values(const std::vector< Point< 2 > > &points, std::vector< std::vector< double > > &values) const override
virtual void vector_gradients(const std::vector< Point< 2 > > &points, std::vector< std::vector< Tensor< 1, 2 > > > &gradients) const override
virtual void vector_laplacians(const std::vector< Point< 2 > > &points, std::vector< std::vector< double > > &values) const override
double lambda() const
The value of lambda.
PoisseuilleFlow(const double r, const double Re)
virtual void vector_gradients(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const override
virtual void vector_values(const std::vector< Point< dim > > &points, std::vector< std::vector< double > > &values) const override
virtual void vector_laplacians(const std::vector< Point< dim > > &points, std::vector< std::vector< double > > &values) const override
virtual void vector_values(const std::vector< Point< dim > > &points, std::vector< std::vector< double > > &values) const override
StokesCosine(const double viscosity=1., const double reaction=0.)
void set_parameters(const double viscosity, const double reaction)
virtual void vector_gradients(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const override
virtual void vector_laplacians(const std::vector< Point< dim > > &points, std::vector< std::vector< double > > &values) const override
virtual void vector_gradients(const std::vector< Point< 2 > > &points, std::vector< std::vector< Tensor< 1, 2 > > > &gradients) const override
virtual void vector_values(const std::vector< Point< 2 > > &points, std::vector< std::vector< double > > &values) const override
double Psi_1(double phi) const
The derivative of Psi()
double Psi(double phi) const
The auxiliary function Psi.
const double lp
Auxiliary variable 1+lambda.
const double lm
Auxiliary variable 1-lambda.
virtual void vector_laplacians(const std::vector< Point< 2 > > &points, std::vector< std::vector< double > > &values) const override
double Psi_3(double phi) const
The 3rd derivative of Psi()
StokesLSingularity()
Constructor setting up some data.
double Psi_4(double phi) const
The 4th derivative of Psi()
double Psi_2(double phi) const
The 2nd derivative of Psi()
const double coslo
Cosine of lambda times omega.
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()
static constexpr double PI
Definition numbers.h:259
STL namespace.
::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 > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)