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)
204 std::vector<std::vector<double>> &
values)
const
206 unsigned int n = points.size();
207 double stretch = 1. / radius;
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;
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;
249 for (
unsigned int d = 0;
d < dim + 1; ++
d)
252 for (
unsigned int k = 0; k < n; ++k)
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)
275 std::vector<std::vector<double>> &
values)
const
277 unsigned int n = points.size();
281 for (
unsigned int d = 0;
d < dim + 1; ++
d)
311 std::vector<std::vector<double>> &
values)
const
313 unsigned int n = points.size();
317 for (
unsigned int d = 0;
d < dim + 1; ++
d)
320 for (
unsigned int k = 0; k < n; ++k)
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;
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;
362 unsigned int n = points.size();
366 for (
unsigned int d = 0;
d < dim + 1; ++
d)
369 for (
unsigned int k = 0; k < n; ++k)
378 const double cx2 = .5 + .5 * c2x;
379 const double cy2 = .5 + .5 * c2y;
395 const double cz2 = .5 + .5 * c2z;
426 std::vector<std::vector<double>> &
values)
const
428 unsigned int n = points.size();
432 for (
unsigned int d = 0;
d < dim + 1; ++
d)
437 vector_values(points,
values);
438 for (
unsigned int d = 0;
d < dim; ++
d)
444 for (
unsigned int d = 0;
d < dim; ++
d)
449 for (
unsigned int k = 0; k < n; ++k)
462 values[0][k] += -viscosity * pi2 * (1. + 2. * c2x) * s2y -
464 values[1][k] += viscosity * pi2 * s2x * (1. + 2. * c2y) -
475 -.5 * viscosity * pi2 * (1. + 2. * c2x) * s2y * s2z -
477 values[1][k] += .5 * viscosity * pi2 * s2x * (1. + 2. * c2y) * s2z -
480 -.5 * viscosity * pi2 * s2x * s2y * (1. + 2. * c2z) -
498 , coslo(
std::
cos(lambda * omega))
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))
568 const double r2 = x * x + y * y;
580 for (
unsigned int d = 0;
d < 3; ++
d)
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))
608 const double r2 = x * x + y * y;
613 const double psi =
Psi(phi);
614 const double psi1 =
Psi_1(phi);
615 const double psi2 =
Psi_2(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)
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)
684 std::vector<std::vector<double>> &
values)
const
686 unsigned int n = points.size();
689 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
692 for (
unsigned int k = 0; k < n; ++k)
695 const double x = p(0);
708 const std::vector<
Point<2>> & points,
711 unsigned int n = points.size();
717 for (
unsigned int i = 0; i < n; ++i)
719 const double x = points[i](0);
720 const double y = points[i](1);
741 std::vector<std::vector<double>> &
values)
const
743 unsigned int n = points.size();
745 for (
unsigned int d = 0;
d < 2 + 1; ++
d)
751 for (
unsigned int k = 0; k < n; ++k)
754 const double x = p(0);
755 const double y = zp * p(1);
757 const double u = 1. - elx *
std::cos(y);
759 const double uy = elx * zp *
std::sin(y);
764 values[0][k] = u * ux + v * uy;
765 values[1][k] = u * vx + v * vy;
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.
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.
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
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)
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)
static constexpr double PI
::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)