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