From 93c6e980fc7fcf58744189a983ba71e25fee6775 Mon Sep 17 00:00:00 2001 From: kanschat Date: Thu, 24 May 2007 00:10:36 +0000 Subject: [PATCH] add Poisseuille flow git-svn-id: https://svn.dealii.org/trunk@14689 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/flow_function.h | 171 ++++++++++++++++ deal.II/base/source/flow_function.cc | 237 ++++++++++++++++++++++ 2 files changed, 408 insertions(+) create mode 100644 deal.II/base/include/base/flow_function.h create mode 100644 deal.II/base/source/flow_function.cc diff --git a/deal.II/base/include/base/flow_function.h b/deal.II/base/include/base/flow_function.h new file mode 100644 index 0000000000..eec8e012da --- /dev/null +++ b/deal.II/base/include/base/flow_function.h @@ -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 +#include +#include + +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 + class FlowFunction : public Function + { + 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 + * p. + */ + 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 >& points, + std::vector >& 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 > &points, + std::vector > > &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 > &points, + std::vector > &values) const = 0; + + virtual void vector_value_list (const std::vector > &points, + std::vector > &values) const; + virtual void vector_gradient_list (const std::vector > &points, + std::vector > > &gradients) const; + /** + * The force term in the + * momentum equation. + */ + virtual void vector_laplacian_list (const std::vector > &points, + std::vector > &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 > aux_values; + + /** + * Auxiliary values for the usual + * Function interface. + */ + mutable std::vector > > aux_gradients; + }; + +/** + * Laminar pipe flow in two and three dimensions. The channel + * stretches along the x-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 + class PoisseuilleFlow : public FlowFunction + { + public: + /** + * Construct an object for the + * given channel radius + * r and the Reynolds + * number Re. + */ + PoisseuilleFlow (const double r, + const double Re); + virtual ~PoisseuilleFlow(); + + 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 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 index 0000000000..166d4f759c --- /dev/null +++ b/deal.II/base/source/flow_function.cc @@ -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 +#include +#include +#include + +#include + +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 + FlowFunction::FlowFunction() + : Function(dim+1), + mean_pressure(0), + aux_values(dim+1), + aux_gradients(dim+1) + {} + + + template + FlowFunction::~FlowFunction() + {} + + + template + void + FlowFunction::pressure_adjustment(double p) + { + mean_pressure = p; + } + + + template + void FlowFunction::vector_value_list ( + const std::vector > &points, + std::vector > &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 + void FlowFunction::vector_gradient_list ( + const std::vector >& points, + std::vector > >& 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 + void FlowFunction::vector_laplacian_list ( + const std::vector >& points, + std::vector >& 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 + unsigned int FlowFunction::memory_consumption () const + { + Assert(false, ExcNotImplemented()); + return 0; + } + + +//----------------------------------------------------------------------// + + template + PoisseuilleFlow::PoisseuilleFlow(const double r, + const double Re) + : + radius(r), Reynolds(Re) + {} + + + template + PoisseuilleFlow::~PoisseuilleFlow() + {} + + + template + void PoisseuilleFlow::vector_values ( + const std::vector >& points, + std::vector >& 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& 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;dmean_pressure; + } + } + + + + template + void PoisseuilleFlow::vector_gradients ( + const std::vector >& points, + std::vector > >& 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& p = points[k]; + // x-velocity + values[0][k][0] = 0.; + for (unsigned int d=1;d + void PoisseuilleFlow::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; + template class FlowFunction<3>; + template class PoisseuilleFlow<2>; + template class PoisseuilleFlow<3>; +} + + + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5