--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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