16 #include <deal.II/base/flow_function.h> 17 #include <deal.II/base/point.h> 18 #include <deal.II/base/tensor.h> 20 #include <deal.II/lac/vector.h> 25 DEAL_II_NAMESPACE_OPEN
35 , aux_gradients(dim + 1)
54 const unsigned int n_points = points.size();
55 Assert(values.size() == n_points,
62 for (
unsigned int d = 0; d < dim + 1; ++d)
66 for (
unsigned int k = 0; k < n_points; ++k)
68 Assert(values[k].size() == dim + 1,
70 for (
unsigned int d = 0; d < dim + 1; ++d)
84 const unsigned int n_points = 1;
85 std::vector<Point<dim>> points(1);
92 for (
unsigned int d = 0; d < dim + 1; ++d)
96 for (
unsigned int d = 0; d < dim + 1; ++d)
104 const unsigned int comp)
const 107 const unsigned int n_points = 1;
108 std::vector<Point<dim>> points(1);
115 for (
unsigned int d = 0; d < dim + 1; ++d)
129 const unsigned int n_points = points.size();
130 Assert(values.size() == n_points,
137 for (
unsigned int d = 0; d < dim + 1; ++d)
141 for (
unsigned int k = 0; k < n_points; ++k)
143 Assert(values[k].size() == dim + 1,
145 for (
unsigned int d = 0; d < dim + 1; ++d)
157 const unsigned int n_points = points.size();
158 Assert(values.size() == n_points,
165 for (
unsigned int d = 0; d < dim + 1; ++d)
169 for (
unsigned int k = 0; k < n_points; ++k)
171 Assert(values[k].size() == dim + 1,
173 for (
unsigned int d = 0; d < dim + 1; ++d)
204 std::vector<std::vector<double>> &values)
const 206 unsigned int n = points.size();
207 double stretch = 1. / radius;
209 Assert(values.size() == dim + 1,
211 for (
unsigned int d = 0; d < dim + 1; ++d)
214 for (
unsigned int k = 0; k < n; ++k)
222 for (
unsigned int d = 1; d < dim; ++d)
223 r2 += p(d) * p(d) * stretch * stretch;
226 values[0][k] = 1. - r2;
228 for (
unsigned int d = 1; d < dim; ++d)
231 values[dim][k] = -2 * (dim - 1) * stretch * stretch * p(0) / Reynolds +
244 unsigned int n = points.size();
245 double stretch = 1. / radius;
247 Assert(values.size() == dim + 1,
249 for (
unsigned int d = 0; d < dim + 1; ++d)
252 for (
unsigned int k = 0; k < n; ++k)
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;
260 for (
unsigned int d = 1; d < dim; ++d)
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.;
275 std::vector<std::vector<double>> &values)
const 277 unsigned int n = points.size();
279 Assert(values.size() == dim + 1,
281 for (
unsigned int d = 0; d < dim + 1; ++d)
284 for (
unsigned int d = 0; d < values.size(); ++d)
285 for (
unsigned int k = 0; k < values[d].size(); ++k)
312 std::vector<std::vector<double>> &values)
const 314 unsigned int n = points.size();
316 Assert(values.size() == dim + 1,
318 for (
unsigned int d = 0; d < dim + 1; ++d)
321 for (
unsigned int k = 0; k < n; ++k)
326 const double cx = cos(x);
327 const double cy = cos(y);
328 const double sx = sin(x);
329 const double sy = sin(y);
333 values[0][k] = cx * cx * cy * sy;
334 values[1][k] = -cx * sx * cy * cy;
340 const double cz = cos(z);
341 const double sz = sin(z);
343 values[0][k] = cx * cx * cy * sy * cz * sz;
344 values[1][k] = cx * sx * cy * cy * cz * sz;
345 values[2][k] = -2. * cx * sx * cy * sy * cz * cz;
346 values[3][k] = cx * sx * cy * sy * cz * sz + this->
mean_pressure;
363 unsigned int n = points.size();
365 Assert(values.size() == dim + 1,
367 for (
unsigned int d = 0; d < dim + 1; ++d)
370 for (
unsigned int k = 0; k < n; ++k)
375 const double c2x = cos(2 * x);
376 const double c2y = cos(2 * y);
377 const double s2x = sin(2 * x);
378 const double s2y = sin(2 * y);
379 const double cx2 = .5 + .5 * c2x;
380 const double cy2 = .5 + .5 * c2y;
394 const double c2z = cos(2 * z);
395 const double s2z = sin(2 * z);
396 const double cz2 = .5 + .5 * c2z;
398 values[0][k][0] = -.125 *
numbers::PI * s2x * s2y * s2z;
399 values[0][k][1] = .25 *
numbers::PI * cx2 * c2y * s2z;
400 values[0][k][2] = .25 *
numbers::PI * cx2 * s2y * c2z;
402 values[1][k][0] = .25 *
numbers::PI * c2x * cy2 * s2z;
403 values[1][k][1] = -.125 *
numbers::PI * s2x * s2y * s2z;
404 values[1][k][2] = .25 *
numbers::PI * s2x * cy2 * c2z;
406 values[2][k][0] = -.5 *
numbers::PI * c2x * s2y * cz2;
407 values[2][k][1] = -.5 *
numbers::PI * s2x * c2y * cz2;
408 values[2][k][2] = .25 *
numbers::PI * s2x * s2y * s2z;
410 values[3][k][0] = .125 *
numbers::PI * c2x * s2y * s2z;
411 values[3][k][1] = .125 *
numbers::PI * s2x * c2y * s2z;
412 values[3][k][2] = .125 *
numbers::PI * s2x * s2y * c2z;
427 std::vector<std::vector<double>> &values)
const 429 unsigned int n = points.size();
431 Assert(values.size() == dim + 1,
433 for (
unsigned int d = 0; d < dim + 1; ++d)
439 for (
unsigned int d = 0; d < dim; ++d)
440 for (
unsigned int k = 0; k < values[d].size(); ++k)
445 for (
unsigned int d = 0; d < dim; ++d)
446 for (
unsigned int k = 0; k < values[d].size(); ++k)
451 for (
unsigned int k = 0; k < n; ++k)
456 const double c2x = cos(2 * x);
457 const double c2y = cos(2 * y);
458 const double s2x = sin(2 * x);
459 const double s2y = sin(2 * y);
464 values[0][k] += -
viscosity * pi2 * (1. + 2. * c2x) * s2y -
465 numbers::PI / 4. * c2x * s2y;
466 values[1][k] +=
viscosity * pi2 * s2x * (1. + 2. * c2y) -
467 numbers::PI / 4. * s2x * c2y;
472 const double z = numbers::PI * p(2);
473 const double c2z = cos(2 * z);
474 const double s2z = sin(2 * z);
477 -.5 *
viscosity * pi2 * (1. + 2. * c2x) * s2y * s2z -
478 numbers::PI / 8. * c2x * s2y * s2z;
479 values[1][k] += .5 *
viscosity * pi2 * s2x * (1. + 2. * c2y) * s2z -
480 numbers::PI / 8. * s2x * c2y * s2z;
482 -.5 *
viscosity * pi2 * s2x * s2y * (1. + 2. * c2z) -
483 numbers::PI / 8. * s2x * s2y * c2z;
500 , coslo(cos(lambda * omega))
509 return coslo * (sin(
lp * phi) /
lp - sin(
lm * phi) /
lm) - cos(
lp * phi) +
517 return coslo * (cos(
lp * phi) - cos(
lm * phi)) +
lp * sin(
lp * phi) -
549 StokesLSingularity::vector_values(
550 const std::vector<
Point<2>> & points,
551 std::vector<std::vector<double>> &values)
const 553 unsigned int n = points.size();
556 for (
unsigned int d = 0; d < 2 + 1; ++d)
559 for (
unsigned int k = 0; k < n; ++k)
562 const double x = p(0);
563 const double y = p(1);
565 if ((x < 0) || (y < 0))
567 const double phi = std::atan2(y, -x) +
numbers::PI;
568 const double r2 = x * x + y * y;
569 const double rl = pow(r2,
lambda / 2.);
570 const double rl1 = pow(r2,
lambda / 2. - .5);
572 rl * (
lp * sin(phi) *
Psi(phi) + cos(phi) *
Psi_1(phi));
574 rl * (
lp * cos(phi) *
Psi(phi) - sin(phi) *
Psi_1(phi));
580 for (
unsigned int d = 0; d < 3; ++d)
589 StokesLSingularity::vector_gradients(
590 const std::vector<
Point<2>> & points,
593 unsigned int n = points.size();
596 for (
unsigned int d = 0; d < 2 + 1; ++d)
599 for (
unsigned int k = 0; k < n; ++k)
602 const double x = p(0);
603 const double y = p(1);
605 if ((x < 0) || (y < 0))
607 const double phi = std::atan2(y, -x) +
numbers::PI;
608 const double r2 = x * x + y * y;
609 const double r = sqrt(r2);
610 const double rl = pow(r2,
lambda / 2.);
611 const double rl1 = pow(r2,
lambda / 2. - .5);
612 const double rl2 = 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 = cos(phi);
617 const double sinp = sin(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);
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);
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;
640 for (
unsigned int d = 0; d < 3; ++d)
649 StokesLSingularity::vector_laplacians(
650 const std::vector<
Point<2>> & points,
651 std::vector<std::vector<double>> &values)
const 653 unsigned int n = points.size();
656 for (
unsigned int d = 0; d < 2 + 1; ++d)
659 for (
unsigned int d = 0; d < values.size(); ++d)
660 for (
unsigned int k = 0; k < values[d].size(); ++k)
671 long double r2 = Reynolds / 2.;
673 long double l = -b / (r2 + std::sqrt(r2 * r2 + b));
678 p_average = 1 / (8 * l) * (std::exp(3. * l) - std::exp(-l));
684 Kovasznay::vector_values(
const std::vector<
Point<2>> & points,
685 std::vector<std::vector<double>> &values)
const 687 unsigned int n = points.size();
690 for (
unsigned int d = 0; d < 2 + 1; ++d)
693 for (
unsigned int k = 0; k < n; ++k)
696 const double x = p(0);
698 const double elx = std::exp(lbda * x);
700 values[0][k] = 1. - elx * cos(y);
701 values[1][k] = .5 /
numbers::PI * lbda * elx * sin(y);
702 values[2][k] = -.5 * elx * elx + p_average + this->
mean_pressure;
708 Kovasznay::vector_gradients(
709 const std::vector<
Point<2>> & points,
710 std::vector<std::vector<
Tensor<1, 2>>> &gradients)
const 712 unsigned int n = points.size();
715 Assert(gradients[0].size() == n,
718 for (
unsigned int i = 0; i < n; ++i)
720 const double x = points[i](0);
721 const double y = points[i](1);
723 const double elx = std::exp(lbda * x);
728 gradients[0][i][0] = -lbda * elx * cy;
730 gradients[1][i][0] = lbda * lbda / (2 *
numbers::PI) * elx * sy;
731 gradients[1][i][1] = lbda * elx * cy;
733 gradients[2][i][0] = -lbda * elx * elx;
734 gradients[2][i][1] = 0.;
741 Kovasznay::vector_laplacians(
const std::vector<
Point<2>> & points,
742 std::vector<std::vector<double>> &values)
const 744 unsigned int n = points.size();
746 for (
unsigned int d = 0; d < 2 + 1; ++d)
752 for (
unsigned int k = 0; k < n; ++k)
755 const double x = p(0);
756 const double y = zp * p(1);
757 const double elx = std::exp(lbda * x);
758 const double u = 1. - elx * cos(y);
759 const double ux = -lbda * elx * cos(y);
760 const double uy = elx * zp * sin(y);
761 const double v = lbda / zp * elx * sin(y);
762 const double vx = lbda * lbda / zp * elx * sin(y);
763 const double vy = zp * lbda / zp * elx * cos(y);
765 values[0][k] = u * ux + v * uy;
766 values[1][k] = u * vx + v * vy;
772 for (
unsigned int d = 0; d < values.size(); ++d)
773 for (
unsigned int k = 0; k < values[d].size(); ++k)
796 DEAL_II_NAMESPACE_CLOSE
std::vector< std::vector< double > > aux_values
virtual void vector_laplacians(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
const double lm
Auxiliary variable 1-lambda.
virtual void vector_laplacians(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
double lambda() const
The value of lambda.
virtual void vector_gradients(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
const double lp
Auxiliary variable 1+lambda.
virtual void vector_laplacians(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const =0
virtual void vector_value_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
static const double lambda
StokesCosine(const double viscosity=1., const double reaction=0.)
virtual void vector_values(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override=0
virtual void vector_gradients(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override=0
void set_parameters(const double viscosity, const double reaction)
double viscosity
The viscosity.
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
double Psi_3(double phi) const
The 3rd derivative of Psi()
static::ExceptionBase & ExcMessage(std::string arg1)
virtual void vector_values(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double reaction
The reaction parameter.
std::vector< std::vector< Tensor< 1, dim > > > aux_gradients
PoisseuilleFlow(const double r, const double Re)
virtual void vector_value(const Point< dim > &points, Vector< double > &value) const override
virtual void vector_values(const std::vector< Point< dim >> &points, std::vector< std::vector< double >> &values) const override
double Psi_2(double phi) const
The 2nd derivative of Psi()
virtual void vector_laplacian_list(const std::vector< Point< dim >> &points, std::vector< Vector< double >> &values) const override
const double coslo
Cosine of lambda times omega.
virtual double value(const Point< dim > &points, const unsigned int component) const override
StokesLSingularity()
Constructor setting up some data.
static::ExceptionBase & ExcNotImplemented()
double Psi(double phi) const
The auxiliary function Psi.
virtual void vector_gradients(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
Kovasznay(const double Re, bool Stokes=false)
void pressure_adjustment(double p)
double Psi_1(double phi) const
The derivative of Psi()
double Psi_4(double phi) const
The 4th derivative of Psi()