l-shaped singularity for Stokes

git-svn-id: https://svn.dealii.org/trunk@14760 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 eec8e012da9d031a1b4bf49b2883858ae3a84c91..c22c8204e3e320152eb2f71994e4f86bea7ebcee 100644 (file)
--- a/deal.II/base/include/base/flow_function.h
+++ b/deal.II/base/include/base/flow_function.h
@@ -163,7 +163,41 @@ namespace Functions
       const double Reynolds;
   };
   
+/**
+ * The solution to Stokes' equations on an L-shaped domain.
+ *
+ * Taken from Houston, Schötzau, Wihler, proceeding ENUMATH 2003.
+ *
+ * @author Guido Kanschat, 2007
+ */
+  class StokesLSingularity : public FlowFunction<2>
+  {
+    public:
+                                      /// Constructor setting upsome data.
+      StokesLSingularity();
       
+      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:
+                                      /// The auxiliary function Psi.
+      double Psi(double phi) const;
+                                      /// The derivative of #Psi
+      double Psi_1(double phi) const;
+                                      /// The 3rd derivative of #Psi
+      double Psi_3(double phi) const;
+                                      /// The angle of the reentrant corner
+      const double omega;
+                                      /// The exponent of the radius
+      static const double lambda;
+                                      /// Cosine of #lambda times #omega
+      const double coslo;
+  };
+  
 }
 
 DEAL_II_NAMESPACE_CLOSE

diff --git a/deal.II/base/source/flow_function.cc b/deal.II/base/source/flow_function.cc
index 166d4f759c0cf91709829409ea815a74a3787786..57805c40c480bd5ebef8194bc1b064a7233fa9da 100644 (file)
--- a/deal.II/base/source/flow_function.cc
+++ b/deal.II/base/source/flow_function.cc
@@ -226,6 +226,111 @@ namespace Functions
        values[d][k] = 0.;
   }
   
+//----------------------------------------------------------------------//
+
+  const double StokesLSingularity::lambda = 0.54448373678246;
+  
+  StokesLSingularity::StokesLSingularity()
+                 :
+                 omega (3./2.*deal_II_numbers::PI),
+                 coslo (std::cos(lambda*omega))
+  {}
+  
+  
+  inline
+  double
+  StokesLSingularity::Psi(double phi) const
+  {
+    return std::sin((1.+lambda) * phi) * coslo / (1.+lambda) - std::cos((1.+lambda) * phi)
+      - std::sin((1.-lambda) * phi) * coslo / (1.-lambda) + std::cos((1.-lambda) * phi);
+  }
+  
+  
+  inline
+  double
+  StokesLSingularity::Psi_1(double phi) const
+  {
+    return std::cos((1.+lambda) * phi) * coslo + (1.+lambda) * std::sin((1.+lambda) * phi)
+      - std::cos((1.-lambda) * phi) * coslo - (1.-lambda) * std::sin((1.-lambda) * phi);
+  }
+  
+  
+  inline
+  double
+  StokesLSingularity::Psi_3(double phi) const
+  {
+    return - (1.+lambda) * (1.+lambda)
+      * (std::cos((1.+lambda) * phi) * coslo + (1.+lambda) * std::sin((1.+lambda) * phi))
+      - (1.-lambda) * (1.-lambda) *
+      (- std::cos((1.-lambda) * phi) * coslo - (1.-lambda) * std::sin((1.-lambda) * phi));
+  }
+  
+  
+  void StokesLSingularity::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 = p(1);
+
+       if ((x<0) || (y<0))
+         {
+           const double phi = std::atan2(y,-x)+M_PI;
+           const double r2 = x*x+y*y;
+           values[0][k] = std::pow(r2,lambda/2.)
+                          * ((1.+lambda) * std::sin(phi) * Psi(phi)
+                             + std::cos(phi) * Psi_1(phi));
+           values[1][k] = std::pow(r2,lambda/2.)
+                          * (std::sin(phi) * Psi_1(phi)
+                             -(1.+lambda) * std::cos(phi) * Psi(phi));
+           values[2][k] = -std::pow(r2,lambda/2.-.5)
+                          * ((1.+lambda) * (1.+lambda) * Psi_1(phi) + Psi_3(phi))
+                          / (1.-lambda);
+         }
+       else
+         {
+           for (unsigned int d=0;d<3;++d)
+             values[d][k] = 0.;
+         }
+      }
+  }
+  
+
+  
+  void StokesLSingularity::vector_gradients (
+    const std::vector<Point<2> >&,
+    std::vector<std::vector<Tensor<1,2> > >&) const
+  {
+    Assert(false, ExcNotImplemented());
+  }
+  
+
+  
+  void StokesLSingularity::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));
+    
+    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>;