From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
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<Point<dim> > &points,
 				      std::vector<std::vector<double> >   &values) const = 0;
 
+      virtual void vector_value (const Point<dim>& points, Vector<double>& value) const;
+      virtual double value (const Point<dim>& points, const unsigned int component) const;
       virtual void vector_value_list (const std::vector<Point<dim> > &points,
 				      std::vector<Vector<double> >   &values) const;
       virtual void vector_gradient_list (const std::vector<Point<dim> >            &points,
@@ -165,6 +167,39 @@ namespace Functions
       const double Reynolds;
   };
   
+  
+/**
+ * Artificial divergence free function with homogeneous boundary
+ * conditions on the cube [-1,1]<sup>dim</sup>.
+ *
+ * @ingroup functions
+ * @author Guido Kanschat, 2007
+ */
+  template <int dim>
+  class StokesCosine :
+      public FlowFunction<dim>
+  {
+    public:
+				       /**
+					* Constructor setting the
+					* Reynolds number required for
+					* pressure computation.
+					*/
+      StokesCosine (const double Reynolds);
+      virtual ~StokesCosine();
+      
+      virtual void vector_values (const std::vector<Point<dim> >& points,
+				  std::vector<std::vector<double> >& values) const;
+      virtual void vector_gradients (const std::vector<Point<dim> >& points,
+				     std::vector<std::vector<Tensor<1,dim> > >& gradients) const;
+      virtual void vector_laplacians (const std::vector<Point<dim> > &points,
+				      std::vector<std::vector<double> >   &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
+ * <i>x=1/2</i>.
+ *
+ * @ingroup functions
+ * @author Guido Kanschat, 2007
+ */
+  class Kovasznay : public FlowFunction<2>
+  {
+    public:
+				       /**
+					* Construct an object for the
+					* give Reynolds number
+					* <tt>Re</tt>. If the
+					* parameter <tt>Stokes</tt> 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<Point<2> >& points,
+				  std::vector<std::vector<double> >& values) const;
+      virtual void vector_gradients (const std::vector<Point<2> >& points,
+				     std::vector<std::vector<Tensor<1,2> > >& gradients) const;
+      virtual void vector_laplacians (const std::vector<Point<2> > &points,
+				      std::vector<std::vector<double> >   &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<int dim>
+  void FlowFunction<dim>::vector_value (
+    const Point<dim>& point,
+    Vector<double>& value) const
+  {
+    Assert(value.size() == dim+1, ExcDimensionMismatch(value.size(), dim+1));
+    
+    const unsigned int n_points = 1;
+    std::vector<Point<dim> > points(1);
+    points[0] = point;
+    
+    for (unsigned int d=0;d<dim+1;++d)
+      aux_values[d].resize(n_points);
+    vector_values(points, aux_values);
+    
+    for (unsigned int d=0;d<dim+1;++d)
+      value(d) = aux_values[d][0];
+  }
+  
+  
+  template<int dim>
+  double FlowFunction<dim>::value (
+    const Point<dim>& point,
+    const unsigned int comp) const
+  {
+    Assert(comp < dim+1, ExcIndexRange(comp, 0, dim+1));
+    const unsigned int n_points = 1;
+    std::vector<Point<dim> > points(1);
+    points[0] = point;
+    
+    for (unsigned int d=0;d<dim+1;++d)
+      aux_values[d].resize(n_points);
+    vector_values(points, aux_values);
+    
+    return aux_values[comp][0];
+  }
+  
+  
   template<int dim>
   void FlowFunction<dim>::vector_gradient_list (
     const std::vector<Point<dim> >& points,
@@ -198,7 +236,7 @@ namespace Functions
 					 // x-velocity
 	values[0][k][0] = 0.;
 	for (unsigned int d=1;d<dim;++d)
-	values[0][k][d] = -2.*p(d)*stretch;
+	values[0][k][d] = -2.*p(d)*stretch*stretch;
 					 // other velocities
 	for (unsigned int d=1;d<dim;++d)
 	  values[d][k] = 0.;	
@@ -226,6 +264,166 @@ namespace Functions
 	values[d][k] = 0.;
   }
   
+//----------------------------------------------------------------------//
+
+  template<int dim>
+  StokesCosine<dim>::StokesCosine(const double Re)
+		  :
+		  Reynolds(Re)
+  {}
+
+  
+  template<int dim>
+  StokesCosine<dim>::~StokesCosine()
+  {}
+
+
+  template<int dim>
+  void StokesCosine<dim>::vector_values (
+    const std::vector<Point<dim> >& points,
+    std::vector<std::vector<double> >& values) const
+  {
+    unsigned int n = points.size();
+    
+    Assert(values.size() == dim+1, ExcDimensionMismatch(values.size(), dim+1));
+    for (unsigned int d=0;d<dim+1;++d)
+      Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
+    
+    for (unsigned int k=0;k<n;++k)
+      {
+	const Point<dim>& 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<int dim>
+  void StokesCosine<dim>::vector_gradients (
+    const std::vector<Point<dim> >& points,
+    std::vector<std::vector<Tensor<1,dim> > >& values) const
+  {
+    unsigned int n = points.size();
+    
+    Assert(values.size() == dim+1, ExcDimensionMismatch(values.size(), dim+1));
+    for (unsigned int d=0;d<dim+1;++d)
+      Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
+    
+    for (unsigned int k=0;k<n;++k)
+      {
+	const Point<dim>& 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<int dim>
+  void StokesCosine<dim>::vector_laplacians (
+    const std::vector<Point<dim> >& points,
+    std::vector<std::vector<double> >& values) const
+  {
+    unsigned int n = points.size();
+    
+    Assert(values.size() == dim+1, ExcDimensionMismatch(values.size(), dim+1));
+    for (unsigned int d=0;d<dim+1;++d)
+      Assert(values[d].size() == n, ExcDimensionMismatch(values[d].size(), n));
+    
+    for (unsigned int k=0;k<n;++k)
+      {
+	const Point<dim>& 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<Point<2> >& points,
+    std::vector<std::vector<double> >& 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<n;++k)
+      {
+	const Point<2>& 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<Point<2> >& points,
+    std::vector<std::vector<Tensor<1,2> > >& 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<n;++i)
+      {
+	const double x = points[i](0);
+	const double y = points[i](1);
+	
+	const double elx = exp(lambda*x);
+	const double cy = cos(2*M_PI*y);
+	const double sy = sin(2*M_PI*y);
+	
+					 // u
+	gradients[0][i][0] = -lambda*elx*cy;
+	gradients[0][i][1] = 2*M_PI*elx*sy;
+	gradients[1][i][0] = lambda*lambda/(2*M_PI)*elx*sy;
+	gradients[1][i][1] =lambda*elx*cy;
+					 // p
+	gradients[2][i][0] = -lambda*elx*elx;
+	gradients[2][i][1] = 0.;
+    }
+  }
+  
+
+  
+  void Kovasznay::vector_laplacians (
+    const std::vector<Point<2> >& points,
+    std::vector<std::vector<double> >& 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<n;++k)
+	  {
+	    const Point<2>& 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<values.size();++d)
+	  for (unsigned int k=0;k<values[d].size();++k)
+	    values[d][k] = 0.;
+      }
+  }
+  
+  
+  
   
   template class FlowFunction<2>;
   template class FlowFunction<3>;
   template class PoisseuilleFlow<2>;
   template class PoisseuilleFlow<3>;
+  template class StokesCosine<2>;
+  template class StokesCosine<3>;
 }