Reference documentation for deal.II version 9.3.3
\(\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\}}\)
flow_function.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2007 - 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
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 : radius(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 unsigned int n = points.size();
207 double stretch = 1. / radius;
208
209 Assert(values.size() == dim + 1,
210 ExcDimensionMismatch(values.size(), dim + 1));
211 for (unsigned int d = 0; d < dim + 1; ++d)
212 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
213
214 for (unsigned int k = 0; k < n; ++k)
215 {
216 const Point<dim> &p = points[k];
217 // First, compute the
218 // square of the distance to
219 // the x-axis divided by the
220 // radius.
221 double r2 = 0;
222 for (unsigned int d = 1; d < dim; ++d)
223 r2 += p(d) * p(d) * stretch * stretch;
224
225 // x-velocity
226 values[0][k] = 1. - r2;
227 // other velocities
228 for (unsigned int d = 1; d < dim; ++d)
229 values[d][k] = 0.;
230 // pressure
231 values[dim][k] = -2 * (dim - 1) * stretch * stretch * p(0) / Reynolds +
232 this->mean_pressure;
233 }
234 }
235
236
237
238 template <int dim>
239 void
241 const std::vector<Point<dim>> & points,
242 std::vector<std::vector<Tensor<1, dim>>> &values) const
243 {
244 unsigned int n = points.size();
245 double stretch = 1. / radius;
246
247 Assert(values.size() == dim + 1,
248 ExcDimensionMismatch(values.size(), dim + 1));
249 for (unsigned int d = 0; d < dim + 1; ++d)
250 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
251
252 for (unsigned int k = 0; k < n; ++k)
253 {
254 const Point<dim> &p = points[k];
255 // x-velocity
256 values[0][k][0] = 0.;
257 for (unsigned int d = 1; d < dim; ++d)
258 values[0][k][d] = -2. * p(d) * stretch * stretch;
259 // other velocities
260 for (unsigned int d = 1; d < dim; ++d)
261 values[d][k] = 0.;
262 // pressure
263 values[dim][k][0] = -2 * (dim - 1) * stretch * stretch / Reynolds;
264 for (unsigned int d = 1; d < dim; ++d)
265 values[dim][k][d] = 0.;
266 }
267 }
268
269
270
271 template <int dim>
272 void
274 const std::vector<Point<dim>> & points,
275 std::vector<std::vector<double>> &values) const
276 {
277 unsigned int n = points.size();
278 (void)n;
279 Assert(values.size() == dim + 1,
280 ExcDimensionMismatch(values.size(), dim + 1));
281 for (unsigned int d = 0; d < dim + 1; ++d)
282 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
283
284 for (auto &point_values : values)
285 std::fill(point_values.begin(), point_values.end(), 0.);
286 }
287
288 //----------------------------------------------------------------------//
289
290 template <int dim>
291 StokesCosine<dim>::StokesCosine(const double nu, const double r)
292 : viscosity(nu)
293 , reaction(r)
294 {}
295
296
297
298 template <int dim>
299 void
300 StokesCosine<dim>::set_parameters(const double nu, const double r)
301 {
302 viscosity = nu;
303 reaction = r;
304 }
305
306
307 template <int dim>
308 void
310 const std::vector<Point<dim>> & points,
311 std::vector<std::vector<double>> &values) const
312 {
313 unsigned int n = points.size();
314
315 Assert(values.size() == dim + 1,
316 ExcDimensionMismatch(values.size(), dim + 1));
317 for (unsigned int d = 0; d < dim + 1; ++d)
318 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
319
320 for (unsigned int k = 0; k < n; ++k)
321 {
322 const Point<dim> &p = points[k];
323 const double x = numbers::PI / 2. * p(0);
324 const double y = numbers::PI / 2. * p(1);
325 const double cx = std::cos(x);
326 const double cy = std::cos(y);
327 const double sx = std::sin(x);
328 const double sy = std::sin(y);
329
330 if (dim == 2)
331 {
332 values[0][k] = cx * cx * cy * sy;
333 values[1][k] = -cx * sx * cy * cy;
334 values[2][k] = cx * sx * cy * sy + this->mean_pressure;
335 }
336 else if (dim == 3)
337 {
338 const double z = numbers::PI / 2. * p(2);
339 const double cz = std::cos(z);
340 const double sz = std::sin(z);
341
342 values[0][k] = cx * cx * cy * sy * cz * sz;
343 values[1][k] = cx * sx * cy * cy * cz * sz;
344 values[2][k] = -2. * cx * sx * cy * sy * cz * cz;
345 values[3][k] = cx * sx * cy * sy * cz * sz + this->mean_pressure;
346 }
347 else
348 {
349 Assert(false, ExcNotImplemented());
350 }
351 }
352 }
353
354
355
356 template <int dim>
357 void
359 const std::vector<Point<dim>> & points,
360 std::vector<std::vector<Tensor<1, dim>>> &values) const
361 {
362 unsigned int n = points.size();
363
364 Assert(values.size() == dim + 1,
365 ExcDimensionMismatch(values.size(), dim + 1));
366 for (unsigned int d = 0; d < dim + 1; ++d)
367 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
368
369 for (unsigned int k = 0; k < n; ++k)
370 {
371 const Point<dim> &p = points[k];
372 const double x = numbers::PI / 2. * p(0);
373 const double y = numbers::PI / 2. * p(1);
374 const double c2x = std::cos(2 * x);
375 const double c2y = std::cos(2 * y);
376 const double s2x = std::sin(2 * x);
377 const double s2y = std::sin(2 * y);
378 const double cx2 = .5 + .5 * c2x; // cos^2 x
379 const double cy2 = .5 + .5 * c2y; // cos^2 y
380
381 if (dim == 2)
382 {
383 values[0][k][0] = -.25 * numbers::PI * s2x * s2y;
384 values[0][k][1] = .5 * numbers::PI * cx2 * c2y;
385 values[1][k][0] = -.5 * numbers::PI * c2x * cy2;
386 values[1][k][1] = .25 * numbers::PI * s2x * s2y;
387 values[2][k][0] = .25 * numbers::PI * c2x * s2y;
388 values[2][k][1] = .25 * numbers::PI * s2x * c2y;
389 }
390 else if (dim == 3)
391 {
392 const double z = numbers::PI / 2. * p(2);
393 const double c2z = std::cos(2 * z);
394 const double s2z = std::sin(2 * z);
395 const double cz2 = .5 + .5 * c2z; // cos^2 z
396
397 values[0][k][0] = -.125 * numbers::PI * s2x * s2y * s2z;
398 values[0][k][1] = .25 * numbers::PI * cx2 * c2y * s2z;
399 values[0][k][2] = .25 * numbers::PI * cx2 * s2y * c2z;
400
401 values[1][k][0] = .25 * numbers::PI * c2x * cy2 * s2z;
402 values[1][k][1] = -.125 * numbers::PI * s2x * s2y * s2z;
403 values[1][k][2] = .25 * numbers::PI * s2x * cy2 * c2z;
404
405 values[2][k][0] = -.5 * numbers::PI * c2x * s2y * cz2;
406 values[2][k][1] = -.5 * numbers::PI * s2x * c2y * cz2;
407 values[2][k][2] = .25 * numbers::PI * s2x * s2y * s2z;
408
409 values[3][k][0] = .125 * numbers::PI * c2x * s2y * s2z;
410 values[3][k][1] = .125 * numbers::PI * s2x * c2y * s2z;
411 values[3][k][2] = .125 * numbers::PI * s2x * s2y * c2z;
412 }
413 else
414 {
415 Assert(false, ExcNotImplemented());
416 }
417 }
418 }
419
420
421
422 template <int dim>
423 void
425 const std::vector<Point<dim>> & points,
426 std::vector<std::vector<double>> &values) const
427 {
428 unsigned int n = points.size();
429
430 Assert(values.size() == dim + 1,
431 ExcDimensionMismatch(values.size(), dim + 1));
432 for (unsigned int d = 0; d < dim + 1; ++d)
433 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
434
435 if (reaction != 0.)
436 {
437 vector_values(points, values);
438 for (unsigned int d = 0; d < dim; ++d)
439 for (double &point_value : values[d])
440 point_value *= -reaction;
441 }
442 else
443 {
444 for (unsigned int d = 0; d < dim; ++d)
445 std::fill(values[d].begin(), values[d].end(), 0.);
446 }
447
448
449 for (unsigned int k = 0; k < n; ++k)
450 {
451 const Point<dim> &p = points[k];
452 const double x = numbers::PI / 2. * p(0);
453 const double y = numbers::PI / 2. * p(1);
454 const double c2x = std::cos(2 * x);
455 const double c2y = std::cos(2 * y);
456 const double s2x = std::sin(2 * x);
457 const double s2y = std::sin(2 * y);
458 const double pi2 = .25 * numbers::PI * numbers::PI;
459
460 if (dim == 2)
461 {
462 values[0][k] += -viscosity * pi2 * (1. + 2. * c2x) * s2y -
463 numbers::PI / 4. * c2x * s2y;
464 values[1][k] += viscosity * pi2 * s2x * (1. + 2. * c2y) -
465 numbers::PI / 4. * s2x * c2y;
466 values[2][k] = 0.;
467 }
468 else if (dim == 3)
469 {
470 const double z = numbers::PI * p(2);
471 const double c2z = std::cos(2 * z);
472 const double s2z = std::sin(2 * z);
473
474 values[0][k] +=
475 -.5 * viscosity * pi2 * (1. + 2. * c2x) * s2y * s2z -
476 numbers::PI / 8. * c2x * s2y * s2z;
477 values[1][k] += .5 * viscosity * pi2 * s2x * (1. + 2. * c2y) * s2z -
478 numbers::PI / 8. * s2x * c2y * s2z;
479 values[2][k] +=
480 -.5 * viscosity * pi2 * s2x * s2y * (1. + 2. * c2z) -
481 numbers::PI / 8. * s2x * s2y * c2z;
482 values[3][k] = 0.;
483 }
484 else
485 {
486 Assert(false, ExcNotImplemented());
487 }
488 }
489 }
490
491
492 //----------------------------------------------------------------------//
493
494 const double StokesLSingularity::lambda = 0.54448373678246;
495
497 : omega(3. / 2. * numbers::PI)
498 , coslo(std::cos(lambda * omega))
499 , lp(1. + lambda)
500 , lm(1. - lambda)
501 {}
502
503
504 inline double
505 StokesLSingularity::Psi(double phi) const
506 {
507 return coslo * (std::sin(lp * phi) / lp - std::sin(lm * phi) / lm) -
508 std::cos(lp * phi) + std::cos(lm * phi);
509 }
510
511
512 inline double
514 {
515 return coslo * (std::cos(lp * phi) - std::cos(lm * phi)) +
516 lp * std::sin(lp * phi) - lm * std::sin(lm * phi);
517 }
518
519
520 inline double
522 {
523 return coslo * (lm * std::sin(lm * phi) - lp * std::sin(lp * phi)) +
524 lp * lp * std::cos(lp * phi) - lm * lm * std::cos(lm * phi);
525 }
526
527
528 inline double
530 {
531 return coslo *
532 (lm * lm * std::cos(lm * phi) - lp * lp * std::cos(lp * phi)) +
533 lm * lm * lm * std::sin(lm * phi) -
534 lp * lp * lp * std::sin(lp * phi);
535 }
536
537
538 inline double
540 {
541 return coslo * (lp * lp * lp * std::sin(lp * phi) -
542 lm * lm * lm * std::sin(lm * phi)) +
543 lm * lm * lm * lm * std::cos(lm * phi) -
544 lp * lp * lp * lp * std::cos(lp * phi);
545 }
546
547
548 void
550 const std::vector<Point<2>> & points,
551 std::vector<std::vector<double>> &values) const
552 {
553 unsigned int n = points.size();
554
555 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
556 for (unsigned int d = 0; d < 2 + 1; ++d)
557 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
558
559 for (unsigned int k = 0; k < n; ++k)
560 {
561 const Point<2> &p = points[k];
562 const double x = p(0);
563 const double y = p(1);
564
565 if ((x < 0) || (y < 0))
566 {
567 const double phi = std::atan2(y, -x) + numbers::PI;
568 const double r2 = x * x + y * y;
569 const double rl = std::pow(r2, lambda / 2.);
570 const double rl1 = std::pow(r2, lambda / 2. - .5);
571 values[0][k] =
572 rl * (lp * std::sin(phi) * Psi(phi) + std::cos(phi) * Psi_1(phi));
573 values[1][k] =
574 rl * (lp * std::cos(phi) * Psi(phi) - std::sin(phi) * Psi_1(phi));
575 values[2][k] = -rl1 * (lp * lp * Psi_1(phi) + Psi_3(phi)) / lm +
576 this->mean_pressure;
577 }
578 else
579 {
580 for (unsigned int d = 0; d < 3; ++d)
581 values[d][k] = 0.;
582 }
583 }
584 }
585
586
587
588 void
590 const std::vector<Point<2>> & points,
591 std::vector<std::vector<Tensor<1, 2>>> &values) const
592 {
593 unsigned int n = points.size();
594
595 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
596 for (unsigned int d = 0; d < 2 + 1; ++d)
597 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
598
599 for (unsigned int k = 0; k < n; ++k)
600 {
601 const Point<2> &p = points[k];
602 const double x = p(0);
603 const double y = p(1);
604
605 if ((x < 0) || (y < 0))
606 {
607 const double phi = std::atan2(y, -x) + numbers::PI;
608 const double r2 = x * x + y * y;
609 const double r = std::sqrt(r2);
610 const double rl = std::pow(r2, lambda / 2.);
611 const double rl1 = std::pow(r2, lambda / 2. - .5);
612 const double rl2 = std::pow(r2, lambda / 2. - 1.);
613 const double psi = Psi(phi);
614 const double psi1 = Psi_1(phi);
615 const double psi2 = Psi_2(phi);
616 const double cosp = std::cos(phi);
617 const double sinp = std::sin(phi);
618
619 // Derivatives of u with respect to r, phi
620 const double udr = lambda * rl1 * (lp * sinp * psi + cosp * psi1);
621 const double udp = rl * (lp * cosp * psi + lp * sinp * psi1 -
622 sinp * psi1 + cosp * psi2);
623 // Derivatives of v with respect to r, phi
624 const double vdr = lambda * rl1 * (lp * cosp * psi - sinp * psi1);
625 const double vdp = rl * (lp * (cosp * psi1 - sinp * psi) -
626 cosp * psi1 - sinp * psi2);
627 // Derivatives of p with respect to r, phi
628 const double pdr =
629 -(lambda - 1.) * rl2 * (lp * lp * psi1 + Psi_3(phi)) / lm;
630 const double pdp = -rl1 * (lp * lp * psi2 + Psi_4(phi)) / lm;
631 values[0][k][0] = cosp * udr - sinp / r * udp;
632 values[0][k][1] = -sinp * udr - cosp / r * udp;
633 values[1][k][0] = cosp * vdr - sinp / r * vdp;
634 values[1][k][1] = -sinp * vdr - cosp / r * vdp;
635 values[2][k][0] = cosp * pdr - sinp / r * pdp;
636 values[2][k][1] = -sinp * pdr - cosp / r * pdp;
637 }
638 else
639 {
640 for (unsigned int d = 0; d < 3; ++d)
641 values[d][k] = 0.;
642 }
643 }
644 }
645
646
647
648 void
650 const std::vector<Point<2>> & points,
651 std::vector<std::vector<double>> &values) const
652 {
653 unsigned int n = points.size();
654 (void)n;
655 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
656 for (unsigned int d = 0; d < 2 + 1; ++d)
657 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
658
659 for (auto &point_values : values)
660 std::fill(point_values.begin(), point_values.end(), 0.);
661 }
662
663
664 //----------------------------------------------------------------------//
665
666 Kovasznay::Kovasznay(double Re, bool stokes)
667 : Reynolds(Re)
668 , stokes(stokes)
669 {
670 long double r2 = Reynolds / 2.;
671 long double b = 4 * numbers::PI * numbers::PI;
672 long double l = -b / (r2 + std::sqrt(r2 * r2 + b));
673 lbda = l;
674 // mean pressure for a domain
675 // spreading from -.5 to 1.5 in
676 // x-direction
677 p_average = 1 / (8 * l) * (std::exp(3. * l) - std::exp(-l));
678 }
679
680
681
682 void
683 Kovasznay::vector_values(const std::vector<Point<2>> & points,
684 std::vector<std::vector<double>> &values) const
685 {
686 unsigned int n = points.size();
687
688 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
689 for (unsigned int d = 0; d < 2 + 1; ++d)
690 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
691
692 for (unsigned int k = 0; k < n; ++k)
693 {
694 const Point<2> &p = points[k];
695 const double x = p(0);
696 const double y = 2. * numbers::PI * p(1);
697 const double elx = std::exp(lbda * x);
698
699 values[0][k] = 1. - elx * std::cos(y);
700 values[1][k] = .5 / numbers::PI * lbda * elx * std::sin(y);
701 values[2][k] = -.5 * elx * elx + p_average + this->mean_pressure;
702 }
703 }
704
705
706 void
708 const std::vector<Point<2>> & points,
709 std::vector<std::vector<Tensor<1, 2>>> &gradients) const
710 {
711 unsigned int n = points.size();
712
713 Assert(gradients.size() == 3, ExcDimensionMismatch(gradients.size(), 3));
714 Assert(gradients[0].size() == n,
715 ExcDimensionMismatch(gradients[0].size(), n));
716
717 for (unsigned int i = 0; i < n; ++i)
718 {
719 const double x = points[i](0);
720 const double y = points[i](1);
721
722 const double elx = std::exp(lbda * x);
723 const double cy = std::cos(2 * numbers::PI * y);
724 const double sy = std::sin(2 * numbers::PI * y);
725
726 // u
727 gradients[0][i][0] = -lbda * elx * cy;
728 gradients[0][i][1] = 2. * numbers::PI * elx * sy;
729 gradients[1][i][0] = lbda * lbda / (2 * numbers::PI) * elx * sy;
730 gradients[1][i][1] = lbda * elx * cy;
731 // p
732 gradients[2][i][0] = -lbda * elx * elx;
733 gradients[2][i][1] = 0.;
734 }
735 }
736
737
738
739 void
740 Kovasznay::vector_laplacians(const std::vector<Point<2>> & points,
741 std::vector<std::vector<double>> &values) const
742 {
743 unsigned int n = points.size();
744 Assert(values.size() == 2 + 1, ExcDimensionMismatch(values.size(), 2 + 1));
745 for (unsigned int d = 0; d < 2 + 1; ++d)
746 Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
747
748 if (stokes)
749 {
750 const double zp = 2. * numbers::PI;
751 for (unsigned int k = 0; k < n; ++k)
752 {
753 const Point<2> &p = points[k];
754 const double x = p(0);
755 const double y = zp * p(1);
756 const double elx = std::exp(lbda * x);
757 const double u = 1. - elx * std::cos(y);
758 const double ux = -lbda * elx * std::cos(y);
759 const double uy = elx * zp * std::sin(y);
760 const double v = lbda / zp * elx * std::sin(y);
761 const double vx = lbda * lbda / zp * elx * std::sin(y);
762 const double vy = zp * lbda / zp * elx * std::cos(y);
763
764 values[0][k] = u * ux + v * uy;
765 values[1][k] = u * vx + v * vy;
766 values[2][k] = 0.;
767 }
768 }
769 else
770 {
771 for (auto &point_values : values)
772 std::fill(point_values.begin(), point_values.end(), 0.);
773 }
774 }
775
776 double
778 {
779 return lbda;
780 }
781
782
783
784 template class FlowFunction<2>;
785 template class FlowFunction<3>;
786 template class PoisseuilleFlow<2>;
787 template class PoisseuilleFlow<3>;
788 template class StokesCosine<2>;
789 template class StokesCosine<3>;
790} // namespace Functions
791
792
793
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.
const double Reynolds
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
static const double lambda
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
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#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 & ExcMessage(std::string arg1)
Expression atan2(const Expression &y, const Expression &x)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::vector< typename FEPointEvaluation< n_components, dim >::value_type > point_values(const Mapping< dim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &vector, const std::vector< Point< spacedim > > &evaluation_points, Utilities::MPI::RemotePointEvaluation< dim, spacedim > &cache, const EvaluationFlags::EvaluationFlags flags=EvaluationFlags::avg)
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
static constexpr double PI
Definition: numbers.h:231
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)