Reference documentation for deal.II version 9.2.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\}}\)
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 
16 #include <deal.II/base/flow_function.h>
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 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 namespace Functions
29 {
30  template <int dim>
31  FlowFunction<dim>::FlowFunction()
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
42  FlowFunction<dim>::pressure_adjustment(double p)
43  {
44  mean_pressure = p;
45  }
46 
47 
48  template <int dim>
49  void
50  FlowFunction<dim>::vector_value_list(
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
78  FlowFunction<dim>::vector_value(const Point<dim> &point,
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
103  FlowFunction<dim>::value(const Point<dim> & point,
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
125  FlowFunction<dim>::vector_gradient_list(
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
153  FlowFunction<dim>::vector_laplacian_list(
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
181  FlowFunction<dim>::memory_consumption() const
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
202  PoisseuilleFlow<dim>::vector_values(
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
240  PoisseuilleFlow<dim>::vector_gradients(
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
273  PoisseuilleFlow<dim>::vector_laplacians(
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
309  StokesCosine<dim>::vector_values(
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
358  StokesCosine<dim>::vector_gradients(
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
424  StokesCosine<dim>::vector_laplacians(
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 
496  StokesLSingularity::StokesLSingularity()
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
513  StokesLSingularity::Psi_1(double phi) const
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
521  StokesLSingularity::Psi_2(double phi) const
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
529  StokesLSingularity::Psi_3(double phi) const
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
539  StokesLSingularity::Psi_4(double phi) const
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
549  StokesLSingularity::vector_values(
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
589  StokesLSingularity::vector_gradients(
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
649  StokesLSingularity::vector_laplacians(
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
707  Kovasznay::vector_gradients(
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
777  Kovasznay::lambda() const
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 
794 DEAL_II_NAMESPACE_CLOSE
Functions::FlowFunction< 2 >::mean_pressure
double mean_pressure
Definition: flow_function.h:125
Functions
Definition: flow_function.h:28
Functions::FlowFunction
Definition: flow_function.h:49
Functions::StokesCosine::vector_gradients
virtual void vector_gradients(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
Definition: flow_function.cc:358
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Functions::StokesLSingularity::Psi_4
double Psi_4(double phi) const
The 4th derivative of Psi()
Definition: flow_function.cc:539
Functions::StokesCosine
Definition: flow_function.h:195
Functions::StokesLSingularity::lambda
static const double lambda
Definition: flow_function.h:284
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
Functions::Kovasznay::lbda
double lbda
Definition: flow_function.h:332
VectorTools::point_value
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
Functions::StokesLSingularity::Psi
double Psi(double phi) const
The auxiliary function Psi.
Definition: flow_function.cc:505
Functions::Kovasznay::lambda
double lambda() const
The value of lambda.
Definition: flow_function.cc:777
Functions::FlowFunction::FlowFunction
FlowFunction()
Definition: flow_function.cc:31
Functions::PoisseuilleFlow::Reynolds
const double Reynolds
Definition: flow_function.h:177
Functions::StokesLSingularity::Psi_3
double Psi_3(double phi) const
The 3rd derivative of Psi()
Definition: flow_function.cc:529
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Functions::Kovasznay::vector_values
virtual void vector_values(const std::vector< Point< 2 >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:683
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
Functions::StokesCosine::vector_values
virtual void vector_values(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:309
Differentiation::SD::atan2
Expression atan2(const Expression &y, const Expression &x)
Definition: symengine_math.cc:154
tensor.h
Functions::StokesLSingularity::Psi_1
double Psi_1(double phi) const
The derivative of Psi()
Definition: flow_function.cc:513
Functions::PoisseuilleFlow::vector_values
virtual void vector_values(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:202
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
Functions::StokesLSingularity::vector_gradients
virtual void vector_gradients(const std::vector< Point< 2 >> &points, std::vector< std::vector< Tensor< 1, 2 >>> &gradients) const override
Definition: flow_function.cc:589
Tensor< 1, dim >
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:268
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Functions::PoisseuilleFlow::vector_laplacians
virtual void vector_laplacians(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:273
Functions::Kovasznay::vector_laplacians
virtual void vector_laplacians(const std::vector< Point< 2 >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:740
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Functions::StokesLSingularity::coslo
const double coslo
Cosine of lambda times omega.
Definition: flow_function.h:286
numbers
Definition: numbers.h:207
flow_function.h
Functions::PoisseuilleFlow
Definition: flow_function.h:153
value
static const bool value
Definition: dof_tools_constraints.cc:433
Functions::PoisseuilleFlow::vector_gradients
virtual void vector_gradients(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
Definition: flow_function.cc:240
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Functions::StokesLSingularity::lp
const double lp
Auxiliary variable 1+lambda.
Definition: flow_function.h:288
Functions::Kovasznay::Kovasznay
Kovasznay(const double Re, bool Stokes=false)
Definition: flow_function.cc:666
Functions::Kovasznay::stokes
const bool stokes
Definition: flow_function.h:334
Functions::StokesCosine::set_parameters
void set_parameters(const double viscosity, const double reaction)
Definition: flow_function.cc:300
vector.h
Functions::StokesLSingularity::StokesLSingularity
StokesLSingularity()
Constructor setting up some data.
Definition: flow_function.cc:496
Functions::Kovasznay::Reynolds
const double Reynolds
Definition: flow_function.h:331
Functions::StokesCosine::StokesCosine
StokesCosine(const double viscosity=1., const double reaction=0.)
Definition: flow_function.cc:291
Functions::StokesCosine::vector_laplacians
virtual void vector_laplacians(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:424
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
Function
Definition: function.h:151
Functions::StokesLSingularity::Psi_2
double Psi_2(double phi) const
The 2nd derivative of Psi()
Definition: flow_function.cc:521
Functions::Kovasznay::p_average
double p_average
Definition: flow_function.h:333
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Functions::StokesLSingularity::lm
const double lm
Auxiliary variable 1-lambda.
Definition: flow_function.h:290
Functions::Kovasznay::vector_gradients
virtual void vector_gradients(const std::vector< Point< 2 >> &points, std::vector< std::vector< Tensor< 1, 2 >>> &gradients) const override
Definition: flow_function.cc:707
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector< double >
Functions::StokesLSingularity::vector_values
virtual void vector_values(const std::vector< Point< 2 >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:549
Functions::StokesLSingularity::vector_laplacians
virtual void vector_laplacians(const std::vector< Point< 2 >> &points, std::vector< std::vector< double >> &values) const override
Definition: flow_function.cc:649
point.h
Generated by   doxygen 1.8.17