From: kanschat Date: Tue, 24 Jul 2007 09:06:05 +0000 (+0000) Subject: new StokesCosine X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4d7e64091fa457e7f04b05a8f81738854542db97;p=dealii-svn.git new StokesCosine git-svn-id: https://svn.dealii.org/trunk@14868 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/flow_function.h b/deal.II/base/include/base/flow_function.h index 29ae0e34f4..2bdde879b2 100644 --- a/deal.II/base/include/base/flow_function.h +++ b/deal.II/base/include/base/flow_function.h @@ -95,6 +95,8 @@ namespace Functions virtual void vector_laplacians (const std::vector > &points, std::vector > &values) const = 0; + virtual void vector_value (const Point& points, Vector& value) const; + virtual double value (const Point& points, const unsigned int component) const; virtual void vector_value_list (const std::vector > &points, std::vector > &values) const; virtual void vector_gradient_list (const std::vector > &points, @@ -165,6 +167,39 @@ namespace Functions const double Reynolds; }; + +/** + * Artificial divergence free function with homogeneous boundary + * conditions on the cube [-1,1]dim. + * + * @ingroup functions + * @author Guido Kanschat, 2007 + */ + template + class StokesCosine : + public FlowFunction + { + public: + /** + * Constructor setting the + * Reynolds number required for + * pressure computation. + */ + StokesCosine (const double Reynolds); + virtual ~StokesCosine(); + + virtual void vector_values (const std::vector >& points, + std::vector >& values) const; + virtual void vector_gradients (const std::vector >& points, + std::vector > >& gradients) const; + virtual void vector_laplacians (const std::vector > &points, + std::vector > &values) const; + + private: + const double Reynolds; + }; + + /** * The solution to Stokes' equations on an L-shaped domain. * @@ -201,6 +236,49 @@ namespace Functions const double coslo; }; +/** + * Flow solution in 2D by Kovasznay (1947). + * + * This function is valid on the half plane right of the line + * x=1/2. + * + * @ingroup functions + * @author Guido Kanschat, 2007 + */ + class Kovasznay : public FlowFunction<2> + { + public: + /** + * Construct an object for the + * give Reynolds number + * Re. If the + * parameter Stokes is + * true, the right hand side of + * the momentum equation + * returned by + * vector_laplacians() contains + * the nonlinearity, such that + * the Kovasznay solution can + * be obtained as the solution + * to a Stokes problem. + */ + Kovasznay (const double Re, bool Stokes = false); + virtual ~Kovasznay(); + + virtual void vector_values (const std::vector >& points, + std::vector >& values) const; + virtual void vector_gradients (const std::vector >& points, + std::vector > >& gradients) const; + virtual void vector_laplacians (const std::vector > &points, + std::vector > &values) const; + + private: + const double Reynolds; + double lambda; + double p_average; + const bool stokes; + }; + } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/base/source/flow_function.cc b/deal.II/base/source/flow_function.cc index 4287941d2c..00c991e6ad 100644 --- a/deal.II/base/source/flow_function.cc +++ b/deal.II/base/source/flow_function.cc @@ -80,6 +80,44 @@ namespace Functions } + template + void FlowFunction::vector_value ( + const Point& point, + Vector& value) const + { + Assert(value.size() == dim+1, ExcDimensionMismatch(value.size(), dim+1)); + + const unsigned int n_points = 1; + std::vector > points(1); + points[0] = point; + + for (unsigned int d=0;d + double FlowFunction::value ( + const Point& point, + const unsigned int comp) const + { + Assert(comp < dim+1, ExcIndexRange(comp, 0, dim+1)); + const unsigned int n_points = 1; + std::vector > points(1); + points[0] = point; + + for (unsigned int d=0;d void FlowFunction::vector_gradient_list ( const std::vector >& points, @@ -198,7 +236,7 @@ namespace Functions // x-velocity values[0][k][0] = 0.; for (unsigned int d=1;d + StokesCosine::StokesCosine(const double Re) + : + Reynolds(Re) + {} + + + template + StokesCosine::~StokesCosine() + {} + + + template + void StokesCosine::vector_values ( + const std::vector >& points, + std::vector >& values) const + { + unsigned int n = points.size(); + + Assert(values.size() == dim+1, ExcDimensionMismatch(values.size(), dim+1)); + for (unsigned int d=0;d& p = points[k]; + const double x = deal_II_numbers::PI * p(0); + const double y = deal_II_numbers::PI * p(1); + const double cx = std::cos(x); + const double cy = std::cos(y); + const double sx = std::sin(x); + const double sy = std::sin(y); + + if (dim==2) + { + values[0][k] = cx*cx*cy*sy; + values[1][k] = -cx*sx*cy*cy; + values[2][k] = cx*sx*cy*sy; + } + else if (dim==3) + { + const double z = deal_II_numbers::PI * p(2); + const double cz = std::cos(z); + const double sz = std::sin(z); + + values[0][k] = cx*cx*cy*sy*cz*sz; + values[1][k] = cx*sx*cy*cy*cz*sz; + values[2][k] = -2.*cx*sx*cy*sy*cz*cz; + values[3][k] = 0.; + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + + + + template + void StokesCosine::vector_gradients ( + const std::vector >& points, + std::vector > >& values) const + { + unsigned int n = points.size(); + + Assert(values.size() == dim+1, ExcDimensionMismatch(values.size(), dim+1)); + for (unsigned int d=0;d& p = points[k]; + const double x = deal_II_numbers::PI * p(0); + const double y = deal_II_numbers::PI * p(1); + const double cx = std::cos(x); + const double cy = std::cos(y); + const double sx = std::sin(x); + const double sy = std::sin(y); + + if (dim==2) + { + values[0][k][0] = -2.*deal_II_numbers::PI * cx*sx*cy*sy; + values[0][k][1] = deal_II_numbers::PI * cx*cx*(cy*cy-sy*sy); + values[1][k][0] = deal_II_numbers::PI * (sx*sx-cx*cx)*cy*cy; + values[1][k][1] = 2.*deal_II_numbers::PI * cx*sx*cy*sy; + values[2][k][0] = deal_II_numbers::PI * (cx*cx-sx*sx)*cy*sy; + values[2][k][1] = deal_II_numbers::PI * cx*sx*(cy*cy-sy*sy); + } + else if (dim==3) + { + const double z = deal_II_numbers::PI * p(2); + const double cz = std::cos(z); + const double sz = std::sin(z); + + values[0][k][0] = -2.*deal_II_numbers::PI * cx*cx*cy*sy*cz*sz; + values[1][k][0] = -2.*deal_II_numbers::PI * cx*sx*cy*cy*cz*sz; + values[2][k][0] = -2.*deal_II_numbers::PI * -2.*cx*sx*cy*sy*cz*cz; + values[3][k][0] = 0.; + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + + + + template + void StokesCosine::vector_laplacians ( + const std::vector >& points, + std::vector >& values) const + { + unsigned int n = points.size(); + + Assert(values.size() == dim+1, ExcDimensionMismatch(values.size(), dim+1)); + for (unsigned int d=0;d& p = points[k]; + const double x = deal_II_numbers::PI * p(0); + const double y = deal_II_numbers::PI * p(1); + const double cx = std::cos(x); + const double cy = std::cos(y); + const double sx = std::sin(x); + const double sy = std::sin(y); + const double prefix = 2. * deal_II_numbers::PI * deal_II_numbers::PI; + + if (dim==2) + { + values[0][k] = prefix * (cx*cx-sx*sx)*cy*sy + + deal_II_numbers::PI * (cx*cx-sx*sx)*cy*sy; + values[1][k] = - prefix * (cy*cy-sy*sy)*cx*sx + + deal_II_numbers::PI * (cy*cy-sy*sy)*cx*sx; + values[2][k] = 0.; + } + else if (dim==3) + { + const double z = deal_II_numbers::PI * p(2); + const double cz = std::cos(z); + const double sz = std::sin(z); + + values[0][k] = cx*cx*cy*sy*cz*sz; + values[1][k] = cx*sx*cy*cy*cz*sz; + values[2][k] = -2.*cx*sx*cy*sy*cz*cz; + values[3][k] = 0.; + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + + //----------------------------------------------------------------------// const double StokesLSingularity::lambda = 0.54448373678246; @@ -330,11 +528,132 @@ namespace Functions } +//----------------------------------------------------------------------// + + Kovasznay::Kovasznay(double Re, bool stokes) + : + Reynolds(Re), + stokes(stokes) + { + long double r2 = Reynolds/2.; + long double b = 4*M_PI*M_PI; + long double l = -b/(r2+sqrt(r2*r2+b)); + lambda = l; + // mean pressure for a domain + // spreading from -.5 to 1.5 in + // x-direction + p_average = 1/(8*l)*(exp(3.*l)-exp(-l)); + } + + + Kovasznay::~Kovasznay() + {} + + + void Kovasznay::vector_values ( + const std::vector >& points, + std::vector >& values) const + { + unsigned int n = points.size(); + + Assert(values.size() == 2+1, ExcDimensionMismatch(values.size(), 2+1)); + for (unsigned int d=0;d<2+1;++d) + Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n)); + + for (unsigned int k=0;k& p = points[k]; + const double x = p(0); + const double y = 2. * deal_II_numbers::PI * p(1); + const double elx = std::exp(lambda*x); + + values[0][k] = 1. - elx * std::cos(y); + values[1][k] = .5 / deal_II_numbers::PI * lambda * elx * std::sin(y); + values[2][k] = .5 * elx * elx - p_average - this->mean_pressure; + } + } + + + void Kovasznay::vector_gradients ( + const std::vector >& points, + std::vector > >& gradients) const + { + Assert(false, ExcNotImplemented()); + unsigned int n = points.size(); + + Assert (gradients.size() == n, ExcDimensionMismatch(gradients.size(), n)); + Assert (gradients[0].size() >= this->n_components, + ExcDimensionMismatch(gradients[0].size(), this->n_components)); + + for (unsigned int i=0;i >& points, + std::vector >& values) const + { + unsigned int n = points.size(); + Assert(values.size() == 2+1, ExcDimensionMismatch(values.size(), 2+1)); + for (unsigned int d=0;d<2+1;++d) + Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n)); + + if (stokes) + { + const double zp = 2. * deal_II_numbers::PI; + for (unsigned int k=0;k& p = points[k]; + const double x = p(0); + const double y = zp * p(1); + const double elx = std::exp(lambda*x); + const double u = 1. - elx * std::cos(y); + const double ux = -lambda * elx * std::cos(y); + const double uy = elx * zp * std::sin(y); + const double v = lambda/zp * elx * std::sin(y); + const double vx = lambda*lambda/zp * elx * std::sin(y); + const double vy = zp*lambda/zp * elx * std::cos(y); + + values[0][k] = u*ux+v*uy; + values[1][k] = u*vx+v*vy; + values[2][k] = 0.; + } + } + else + { + for (unsigned int d=0;d; template class FlowFunction<3>; template class PoisseuilleFlow<2>; template class PoisseuilleFlow<3>; + template class StokesCosine<2>; + template class StokesCosine<3>; }