From cb2871dc4c79c14c4edaf23a304c0836a491adb5 Mon Sep 17 00:00:00 2001 From: guido Date: Mon, 12 Mar 2001 00:07:28 +0000 Subject: [PATCH] variable directions git-svn-id: https://svn.dealii.org/trunk@4183 0785d39b-7218-0410-832d-ea1e28bc413d --- .../base/include/base/function_derivative.h | 45 +++++++++----- deal.II/base/source/function_derivative.cc | 59 +++++++++++++------ 2 files changed, 70 insertions(+), 34 deletions(-) diff --git a/deal.II/base/include/base/function_derivative.h b/deal.II/base/include/base/function_derivative.h index bffcabc2c7..714a65af85 100644 --- a/deal.II/base/include/base/function_derivative.h +++ b/deal.II/base/include/base/function_derivative.h @@ -49,19 +49,39 @@ class FunctionDerivative : public Function /** * Constructor. Provided are the - * function to compute derivatives - * of and the direction vector of - * the differentiation. + * function to compute + * derivatives of, the direction + * vector of the differentiation + * and the step size @p{h} of the + * difference formula. */ FunctionDerivative (const Function& f, - const Point& direction); + const Point& direction, + const double h = 1.e-6); /** - * Set step size of the difference - * formula. This is set to the - * default in the constructor. + * Constructor. Provided are the + * function to compute + * derivatives of and the + * direction vector of the + * differentiation in each + * quadrature point and the + * difference step size. + * + * This is the constructor for a + * variable velocity field. Most + * probably, a new object of + * @p{FunctionDerivative} has to + * be constructed for each set of + * quadrature points. + * + * The number of quadrature point + * must still be the same, when + * values are accessed. */ - void set_h (double h_new = 1.e-6); + FunctionDerivative (const Function& f, + const vector >& direction, + const double h = 1.e-6); /** * Choose the difference formula. @@ -109,17 +129,12 @@ class FunctionDerivative : public Function unsigned int memory_consumption () const; DeclException0(ExcInvalidFormula); - private: /** * Function for differentiation. */ const Function& f; - /** - * Differentiation vector. - */ - Point direction; - + /** * Step size of the difference formula. */ @@ -133,7 +148,7 @@ class FunctionDerivative : public Function * increment vector for the * formula. */ - Point incr; + vector > incr; }; #endif diff --git a/deal.II/base/source/function_derivative.cc b/deal.II/base/source/function_derivative.cc index be2607c6f4..ada97de275 100644 --- a/deal.II/base/source/function_derivative.cc +++ b/deal.II/base/source/function_derivative.cc @@ -19,24 +19,32 @@ template FunctionDerivative::FunctionDerivative (const Function& f, - const Point& dir) + const Point& dir, + const double h) : Function (f.n_components, f.get_time()), f(f), - direction(dir) + h(h), + incr(1, h*dir) { - set_h(); set_formula(); } template -void -FunctionDerivative::set_h (double h_new) +FunctionDerivative::FunctionDerivative (const Function& f, + const vector >& dir, + const double h) + : + Function (f.n_components, f.get_time()), + f(f), + h(h), + incr(dir.size()) { - h = h_new; - incr = h*direction; + for (unsigned int i=0;i::value (const Point &p, const unsigned int component) const { + Assert (incr.size() == 1, + ExcMessage ("FunctionDerivative was not initialized for constant direection")); + switch (formula) { case Euler: - return (f.value(p+incr, component)-f.value(p-incr, component))/(2*h); + return (f.value(p+incr[0], component)-f.value(p-incr[0], component))/(2*h); case UpwindEuler: - return (f.value(p, component)-f.value(p-incr, component))/h; + return (f.value(p, component)-f.value(p-incr[0], component))/h; case FourthOrder: - return (-f.value(p+2*incr, component) + 8*f.value(p+incr, component) - -8*f.value(p-incr, component) + f.value(p-2*incr, component))/(12*h); + return (-f.value(p+2*incr[0], component) + 8*f.value(p+incr[0], component) + -8*f.value(p-incr[0], component) + f.value(p-2*incr[0], component))/(12*h); default: Assert(false, ExcInvalidFormula()); } @@ -77,6 +88,7 @@ FunctionDerivative::value (const Point &p, using namespace std; #endif +//TODO:Optimize construction of vectors! template void @@ -85,7 +97,11 @@ FunctionDerivative::value_list (const typename std::vector > &po const unsigned int component) const { const unsigned int n = points.size(); - + bool variable_direction = (incr.size() == 1) ? false : true; + if (variable_direction) + Assert (incr.size() == points.size(), + ExcDimensionMismatch(incr.size(), points.size())); + switch (formula) { case Euler: @@ -95,8 +111,9 @@ FunctionDerivative::value_list (const typename std::vector > &po std::vector e2(n); for (unsigned int i=0;i::value_list (const typename std::vector > &po std::vector > p2(n); std::vector e2(n); for (unsigned int i=0;i::value_list (const typename std::vector > &po std::vector e_m(n); for (unsigned int i=0;i