Reference documentation for deal.II version GIT d1e6af55b6 2023-02-07 07:50:02+00:00
\(\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 
28 namespace 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);
114 
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  }
121 
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
510  StokesLSingularity::Psi_1(double phi) const
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
518  StokesLSingularity::Psi_2(double phi) const
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
526  StokesLSingularity::Psi_3(double phi) const
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
536  StokesLSingularity::Psi_4(double phi) const
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 
void pressure_adjustment(double p)
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
virtual void vector_gradient_list(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
virtual std::size_t memory_consumption() const override
virtual void vector_value(const Point< dim > &points, Vector< double > &value) const override
virtual void vector_laplacian_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) 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_laplacians(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
double lambda() const
The value of lambda.
const double Reynolds
PoisseuilleFlow(const double r, const double Re)
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_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
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_laplacians(const std::vector< Point< 2 >> &points, std::vector< std::vector< double >> &values) const override
static const double lambda
double Psi_1(double phi) const
The derivative of Psi()
virtual void vector_values(const std::vector< Point< 2 >> &points, std::vector< std::vector< double >> &values) const override
double Psi(double phi) const
The auxiliary function Psi.
const double lp
Auxiliary variable 1+lambda.
virtual void vector_gradients(const std::vector< Point< 2 >> &points, std::vector< std::vector< Tensor< 1, 2 >>> &gradients) const override
const double lm
Auxiliary variable 1-lambda.
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:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1821
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:189
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
std::vector< typename FEPointEvaluation< n_components, dim, spacedim, typename VectorType::value_type >::value_type > point_values(const Mapping< dim > &mapping, const MeshType< dim, spacedim > &mesh, const VectorType &vector, const std::vector< Point< spacedim >> &evaluation_points, Utilities::MPI::RemotePointEvaluation< dim, spacedim > &cache, const EvaluationFlags::EvaluationFlags flags=EvaluationFlags::avg, const unsigned int first_selected_component=0)
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:248