]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add Poisseuille flow
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 24 May 2007 00:10:36 +0000 (00:10 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 24 May 2007 00:10:36 +0000 (00:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@14689 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/flow_function.h [new file with mode: 0644]
deal.II/base/source/flow_function.cc [new file with mode: 0644]

diff --git a/deal.II/base/include/base/flow_function.h b/deal.II/base/include/base/flow_function.h
new file mode 100644 (file)
index 0000000..eec8e01
--- /dev/null
@@ -0,0 +1,171 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__flow_function_h
+#define __deal2__flow_function_h
+
+
+#include <base/config.h>
+#include <base/function.h>
+#include <base/point.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+/**
+ * Base class for analytic solutions to incompressible flow problems.
+ *
+ * Additional to the Function interface, this function provides for an
+ * offset of the pressure: if the pressure of the computed solution
+ * has an integral mean value different from zero, this value can be
+ * given to pressure_adjustment() in order to compute correct pressure
+ * errors.
+ *
+ * @note Derived classes should implement pressures with integral mean
+ * value zero always.
+ *
+ * @note Thread safety: Some of the functions make use of internal data to compute
+ * values. Therefore, every thread should obtain its own object of
+ * derived classes.
+ *
+ * @author Guido Kanschat, 2007
+ */
+  template <int dim>
+  class FlowFunction : public Function<dim>
+  {
+    public:
+                                      /**
+                                       * Constructor, setting up some
+                                       * internal data structures.
+                                       */
+      FlowFunction();
+
+                                      /**
+                                       * Virtual destructor.
+                                       */
+      virtual ~FlowFunction();
+      
+                                      /**
+                                       * Store an adjustment for the
+                                       * pressure function, such that
+                                       * its mean value is
+                                       * <tt>p</tt>.
+                                       */
+      void pressure_adjustment(double p);
+
+                                      /**
+                                       * Values in a structure more
+                                       * suitable for vector valued
+                                       * functions. The outer vector
+                                       * is indexed by solution
+                                       * component, the inner by
+                                       * quadrature point.
+                                       */
+      virtual void vector_values (const std::vector<Point<dim> >& points,
+                                 std::vector<std::vector<double> >& values) const = 0;
+                                      /**
+                                       * Gradients in a structure more
+                                       * suitable for vector valued
+                                       * functions. The outer vector
+                                       * is indexed by solution
+                                       * component, the inner by
+                                       * quadrature point.
+                                       */
+      virtual void vector_gradients (const std::vector<Point<dim> >            &points,
+                                    std::vector<std::vector<Tensor<1,dim> > > &gradients) const = 0;
+                                      /**
+                                       * Force terms in a structure more
+                                       * suitable for vector valued
+                                       * functions. The outer vector
+                                       * is indexed by solution
+                                       * component, the inner by
+                                       * quadrature point.
+                                       */
+      virtual void vector_laplacians (const std::vector<Point<dim> > &points,
+                                     std::vector<std::vector<double> >   &values) const = 0;
+
+      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,
+                                        std::vector<std::vector<Tensor<1,dim> > > &gradients) const;
+                                      /**
+                                       * The force term in the
+                                       * momentum equation.
+                                       */
+      virtual void vector_laplacian_list (const std::vector<Point<dim> > &points,
+                                         std::vector<Vector<double> >   &values) const;
+      
+      unsigned int memory_consumption () const;
+
+    protected:
+                                      /**
+                                       * Mean value of the pressure
+                                       * to be added by derived
+                                       * classes.
+                                       */
+      double mean_pressure;
+
+    private:  
+                                      /**
+                                       * Auxiliary values for the usual
+                                       * Function interface.
+                                       */
+      mutable std::vector<std::vector<double> > aux_values;
+      
+                                      /**
+                                       * Auxiliary values for the usual
+                                       * Function interface.
+                                       */
+      mutable std::vector<std::vector<Tensor<1,dim> > > aux_gradients;
+  };
+
+/**
+ * Laminar pipe flow in two and three dimensions. The channel
+ * stretches along the <i>x</i>-axis and has radius #radius. The
+ * #Reynolds number is used to scale the pressure properly for a
+ * Navier-Stokes problem.
+ *
+ * @author Guido Kanschat, 2007
+ */
+  template <int dim>
+  class PoisseuilleFlow : public FlowFunction<dim>
+  {
+    public:
+                                      /**
+                                       * Construct an object for the
+                                       * given channel radius
+                                       * <tt>r</tt> and the Reynolds
+                                       * number <tt>Re</tt>.
+                                       */
+      PoisseuilleFlow<dim> (const double r,
+                           const double Re);
+      virtual ~PoisseuilleFlow();
+      
+      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 radius;
+      const double Reynolds;
+  };
+  
+      
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/base/source/flow_function.cc b/deal.II/base/source/flow_function.cc
new file mode 100644 (file)
index 0000000..166d4f7
--- /dev/null
@@ -0,0 +1,237 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <base/tensor.h>
+#include <base/point.h>
+#include <base/flow_function.h>
+#include <lac/vector.h>
+
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+// in strict ANSI C mode, the following constants are not defined by
+// default, so we do it ourselves
+#ifndef M_PI
+#  define      M_PI            3.14159265358979323846
+#endif
+
+#ifndef M_PI_2
+#  define      M_PI_2          1.57079632679489661923
+#endif
+
+
+
+namespace Functions
+{
+  
+  template<int dim>
+  FlowFunction<dim>::FlowFunction()
+                 : Function<dim>(dim+1),
+                   mean_pressure(0),
+                   aux_values(dim+1),
+                   aux_gradients(dim+1)
+  {}
+
+  
+  template<int dim>
+  FlowFunction<dim>::~FlowFunction()
+  {}
+
+  
+  template<int dim>
+  void
+  FlowFunction<dim>::pressure_adjustment(double p)
+  {
+    mean_pressure = p;
+  }
+
+  
+  template<int dim>
+  void FlowFunction<dim>::vector_value_list (
+    const std::vector<Point<dim> > &points,
+    std::vector<Vector<double> >   &values) const
+  {
+    const unsigned int n_points = points.size();
+    Assert(values.size() == n_points, ExcDimensionMismatch(values.size(), n_points));
+    
+    for (unsigned int d=0;d<dim+1;++d)
+      aux_values[d].resize(n_points);
+    vector_values(points, aux_values);
+
+    for (unsigned int k=0;k<n_points;++k)
+      {
+       Assert(values[k].size() == dim+1, ExcDimensionMismatch(values[k].size(), dim+1));
+       for (unsigned int d=0;d<dim+1;++d)
+         values[k](d) = aux_values[d][k];
+      }
+  }
+
+  
+  template<int dim>
+  void FlowFunction<dim>::vector_gradient_list (
+    const std::vector<Point<dim> >& points,
+    std::vector<std::vector<Tensor<1,dim> > >& values) const
+  {
+    const unsigned int n_points = points.size();
+    Assert(values.size() == n_points, ExcDimensionMismatch(values.size(), n_points));
+    
+    for (unsigned int d=0;d<dim+1;++d)
+      aux_gradients[d].resize(n_points);
+    vector_gradients(points, aux_gradients);
+
+    for (unsigned int k=0;k<n_points;++k)
+      {
+       Assert(values[k].size() == dim+1, ExcDimensionMismatch(values[k].size(), dim+1));
+       for (unsigned int d=0;d<dim+1;++d)
+         values[k][d] = aux_gradients[d][k];
+      }
+  }
+  
+  
+  template<int dim>
+  void FlowFunction<dim>::vector_laplacian_list (
+    const std::vector<Point<dim> >& points,
+    std::vector<Vector<double> >& values) const
+  {
+    const unsigned int n_points = points.size();
+    Assert(values.size() == n_points, ExcDimensionMismatch(values.size(), n_points));
+    
+    for (unsigned int d=0;d<dim+1;++d)
+      aux_values[d].resize(n_points);
+    vector_laplacians(points, aux_values);
+
+    for (unsigned int k=0;k<n_points;++k)
+      {
+       Assert(values[k].size() == dim+1, ExcDimensionMismatch(values[k].size(), dim+1));
+       for (unsigned int d=0;d<dim+1;++d)
+         values[k](d) = aux_values[d][k];
+      }
+  }
+
+  
+  template<int dim>
+  unsigned int FlowFunction<dim>::memory_consumption () const
+  {
+    Assert(false, ExcNotImplemented());
+    return 0;
+  }
+  
+
+//----------------------------------------------------------------------//
+
+  template<int dim>
+  PoisseuilleFlow<dim>::PoisseuilleFlow(const double r,
+                                       const double Re)
+                 :
+                 radius(r), Reynolds(Re)
+  {}
+
+  
+  template<int dim>
+  PoisseuilleFlow<dim>::~PoisseuilleFlow()
+  {}
+
+
+  template<int dim>
+  void PoisseuilleFlow<dim>::vector_values (
+    const std::vector<Point<dim> >& points,
+    std::vector<std::vector<double> >& values) const
+  {
+    unsigned int n = points.size();
+    double stretch = 1./radius;
+
+    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];
+                                        // First, compute the
+                                        // square of the distance to
+                                        // the x-axis divided by the
+                                        // radius.
+       double r2 = 0;
+       for (unsigned int d=1;d<dim;++d)
+         r2 += p(d)*p(d)*stretch*stretch;
+
+                                        // x-velocity
+       values[0][k] = 1.-r2;
+                                        // other velocities
+       for (unsigned int d=1;d<dim;++d)
+         values[d][k] = 0.;
+                                        // pressure
+       values[dim][k] = -2*(dim-1)*stretch*stretch*p(0)/Reynolds+this->mean_pressure;
+      }
+  }
+  
+
+  
+  template<int dim>
+  void PoisseuilleFlow<dim>::vector_gradients (
+    const std::vector<Point<dim> >& points,
+    std::vector<std::vector<Tensor<1,dim> > >& values) const
+  {
+    unsigned int n = points.size();
+    double stretch = 1./radius;
+
+    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];
+                                        // x-velocity
+       values[0][k][0] = 0.;
+       for (unsigned int d=1;d<dim;++d)
+       values[0][k][d] = -2.*p(d)*stretch;
+                                        // other velocities
+       for (unsigned int d=1;d<dim;++d)
+         values[d][k] = 0.;    
+                                        // pressure
+       values[dim][k][0] = -2*(dim-1)*stretch*stretch/Reynolds;
+       for (unsigned int d=1;d<dim;++d)
+         values[dim][k][d] = 0.;
+      }
+  }
+  
+
+  
+  template<int dim>
+  void PoisseuilleFlow<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 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>;
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.