35 , aux_gradients(dim + 1)
54 const unsigned int n_points = points.size();
60 std::lock_guard<std::mutex> lock(mutex);
62 for (
unsigned int d = 0;
d < dim + 1; ++
d)
63 aux_values[
d].resize(n_points);
64 vector_values(points, aux_values);
66 for (
unsigned int k = 0; k < n_points; ++k)
70 for (
unsigned int d = 0;
d < dim + 1; ++
d)
81 Assert(value.size() == dim + 1,
84 const unsigned int n_points = 1;
85 std::vector<Point<dim>> points(1);
90 std::lock_guard<std::mutex> lock(mutex);
92 for (
unsigned int d = 0;
d < dim + 1; ++
d)
93 aux_values[
d].resize(n_points);
94 vector_values(points, aux_values);
96 for (
unsigned int d = 0;
d < dim + 1; ++
d)
97 value(
d) = aux_values[
d][0];
104 const unsigned int comp)
const
107 const unsigned int n_points = 1;
108 std::vector<Point<dim>> points(1);
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);
119 return aux_values[comp][0];
129 const unsigned int n_points = points.size();
135 std::lock_guard<std::mutex> lock(mutex);
137 for (
unsigned int d = 0;
d < dim + 1; ++
d)
138 aux_gradients[
d].resize(n_points);
139 vector_gradients(points, aux_gradients);
141 for (
unsigned int k = 0; k < n_points; ++k)
145 for (
unsigned int d = 0;
d < dim + 1; ++
d)
157 const unsigned int n_points = points.size();
163 std::lock_guard<std::mutex> lock(mutex);
165 for (
unsigned int d = 0;
d < dim + 1; ++
d)
166 aux_values[
d].resize(n_points);
167 vector_laplacians(points, aux_values);
169 for (
unsigned int k = 0; k < n_points; ++k)
173 for (
unsigned int d = 0;
d < dim + 1; ++
d)
192 : inv_sqr_radius(1 / r / r)
204 std::vector<std::vector<double>> &
values)
const
206 const unsigned int n = points.size();
210 for (
unsigned int d = 0;
d < dim + 1; ++
d)
213 for (
unsigned int k = 0; k < n; ++k)
219 for (
unsigned int d = 1;
d < dim; ++
d)
221 r2 *= inv_sqr_radius;
226 for (
unsigned int d = 1;
d < dim; ++
d)
229 values[dim][k] = -2 * (dim - 1) * inv_sqr_radius * p(0) / Reynolds +
242 const unsigned int n = points.size();
246 for (
unsigned int d = 0;
d < dim + 1; ++
d)
249 for (
unsigned int k = 0; k < n; ++k)
254 for (
unsigned int d = 1;
d < dim; ++
d)
255 values[0][k][
d] = -2. * p(
d) * inv_sqr_radius;
257 for (
unsigned int d = 1;
d < dim; ++
d)
260 values[dim][k][0] = -2 * (dim - 1) * inv_sqr_radius / Reynolds;
261 for (
unsigned int d = 1;
d < dim; ++
d)
272 std::vector<std::vector<double>> &
values)
const
274 const unsigned int n = points.size();
278 for (
unsigned int d = 0;
d < dim + 1; ++
d)
308 std::vector<std::vector<double>> &
values)
const
310 unsigned int n = points.size();
314 for (
unsigned int d = 0;
d < dim + 1; ++
d)
317 for (
unsigned int k = 0; k < n; ++k)
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);
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;
336 const double cz = std::cos(z);
337 const double sz = std::sin(z);
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;
359 unsigned int n = points.size();
363 for (
unsigned int d = 0;
d < dim + 1; ++
d)
366 for (
unsigned int k = 0; k < n; ++k)
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;
376 const double cy2 = .5 + .5 * c2y;
390 const double c2z = std::cos(2 * z);
391 const double s2z = std::sin(2 * z);
392 const double cz2 = .5 + .5 * c2z;
423 std::vector<std::vector<double>> &
values)
const
425 unsigned int n = points.size();
429 for (
unsigned int d = 0;
d < dim + 1; ++
d)
434 vector_values(points,
values);
435 for (
unsigned int d = 0;
d < dim; ++
d)
441 for (
unsigned int d = 0;
d < dim; ++
d)
446 for (
unsigned int k = 0; k < n; ++k)
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);
459 values[0][k] += -viscosity * pi2 * (1. + 2. * c2x) * s2y -
461 values[1][k] += viscosity * pi2 * s2x * (1. + 2. * c2y) -
468 const double c2z = std::cos(2 * z);
469 const double s2z = std::sin(2 * z);
472 -.5 * viscosity * pi2 * (1. + 2. * c2x) * s2y * s2z -
474 values[1][k] += .5 * viscosity * pi2 * s2x * (1. + 2. * c2y) * s2z -
477 -.5 * viscosity * pi2 * s2x * s2y * (1. + 2. * c2z) -
495 , coslo(std::cos(lambda * omega))
504 return coslo * (std::sin(
lp * phi) /
lp - std::sin(
lm * phi) /
lm) -
505 std::cos(
lp * phi) + std::cos(
lm * phi);
512 return coslo * (std::cos(
lp * phi) - std::cos(
lm * phi)) +
513 lp * std::sin(
lp * phi) -
lm * std::sin(
lm * phi);
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);
529 (
lm *
lm * std::cos(
lm * phi) -
lp *
lp * std::cos(
lp * phi)) +
539 lm *
lm *
lm * std::sin(
lm * phi)) +
547 const std::vector<
Point<2>> &points,
548 std::vector<std::vector<double>> &
values)
const
550 unsigned int n = points.size();
553 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
556 for (
unsigned int k = 0; k < n; ++k)
559 const double x = p(0);
560 const double y = p(1);
562 if ((x < 0) || (y < 0))
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);
569 rl * (
lp * std::sin(phi) *
Psi(phi) + std::cos(phi) *
Psi_1(phi));
571 rl * (
lp * std::cos(phi) *
Psi(phi) - std::sin(phi) *
Psi_1(phi));
577 for (
unsigned int d = 0;
d < 3; ++
d)
587 const std::vector<
Point<2>> &points,
590 unsigned int n = points.size();
593 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
596 for (
unsigned int k = 0; k < n; ++k)
599 const double x = p(0);
600 const double y = p(1);
602 if ((x < 0) || (y < 0))
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);
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);
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);
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;
637 for (
unsigned int d = 0;
d < 3; ++
d)
647 const std::vector<
Point<2>> &points,
648 std::vector<std::vector<double>> &
values)
const
650 unsigned int n = points.size();
653 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
669 long double l = -
b / (r2 + std::sqrt(r2 * r2 +
b));
674 p_average = 1 / (8 *
l) * (std::exp(3. *
l) - std::exp(-
l));
681 std::vector<std::vector<double>> &
values)
const
683 unsigned int n = points.size();
686 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
689 for (
unsigned int k = 0; k < n; ++k)
692 const double x = p(0);
694 const double elx = std::exp(
lbda * x);
696 values[0][k] = 1. - elx * std::cos(y);
705 const std::vector<
Point<2>> &points,
708 unsigned int n = points.size();
714 for (
unsigned int i = 0; i < n; ++i)
716 const double x = points[i](0);
717 const double y = points[i](1);
719 const double elx = std::exp(
lbda * x);
738 std::vector<std::vector<double>> &
values)
const
740 unsigned int n = points.size();
742 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
748 for (
unsigned int k = 0; k < n; ++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);
761 values[0][k] = u * ux + v * uy;
762 values[1][k] = u * vx + v * vy;
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.
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.
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
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)
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)
static constexpr double PI