From 820de33707c0e75325e519d9d9f4ff22426df70f Mon Sep 17 00:00:00 2001 From: guido Date: Sun, 12 Nov 2000 21:36:25 +0000 Subject: [PATCH] FunctionDerivative: probably buggy git-svn-id: https://svn.dealii.org/trunk@3479 0785d39b-7218-0410-832d-ea1e28bc413d --- .../base/include/base/function_derivative.h | 111 ++++++++++++++++ deal.II/base/source/function_derivative.cc | 118 ++++++++++++++++++ 2 files changed, 229 insertions(+) create mode 100644 deal.II/base/include/base/function_derivative.h create mode 100644 deal.II/base/source/function_derivative.cc diff --git a/deal.II/base/include/base/function_derivative.h b/deal.II/base/include/base/function_derivative.h new file mode 100644 index 0000000000..e080a5ac4c --- /dev/null +++ b/deal.II/base/include/base/function_derivative.h @@ -0,0 +1,111 @@ +//------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 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__function_derivative_h +#define __deal2__function_derivative_h + +#include +#include + +/** + * Names of difference formulas. + */ +enum DifferenceFormula +{ + Euler, + UpwindEuler +}; + + +/** + * Derivative of a function object. The value access functions of + * this class return the directional derivative of a function with + * respect to a direction provided on construction. If @p{b} is the + * vector, the derivative @p{b . grad f} is computed. This derivative + * is evaluated directly, not by computing the gradient of @p{f} and + * its scalar product with @p{b}. + * + * The derivative is computed numerically, using one of the provided + * difference formulas. + * + * @author Guido Kanschat, 2000 + */ +template +class FunctionDerivative : public Function +{ +public: + /** + * Constructor. Provided are the + * function to compute derivatives + * of and the direction vector of + * the differentiation. + */ + FunctionDerivative (const Function& f, + const Point& direction); + + /** + * Set step size of the difference + * formula. This is set to the + * default in the constructor. + */ + void set_h (double h_new = 1.e-6); + + /** + * Choose the difference formula. + * This is set to the default in + * the constructor. + */ + void set_formula (DifferenceFormula formula = Euler); + + /** + * Function value at one point. + */ + virtual double value (const Point &p, + const unsigned int component = 0) const; + + /** + * Function values at multiple points. + */ + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; + + + DeclException0(ExcInvalidFormula); + +private: + /** + * Function for differentiation. + */ + const Function& f; + /** + * Differentiation vector. + */ + Point direction; + + /** + * Step size of the difference formula. + */ + double h; + /** + * Difference formula. + */ + DifferenceFormula formula; + /** + * Helper object. Contains the + * increment vector for the + * formula. + */ + Point incr; +}; + +#endif diff --git a/deal.II/base/source/function_derivative.cc b/deal.II/base/source/function_derivative.cc new file mode 100644 index 0000000000..5866eb9747 --- /dev/null +++ b/deal.II/base/source/function_derivative.cc @@ -0,0 +1,118 @@ +//-------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 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 + +template +FunctionDerivative::FunctionDerivative (const Function& f, + const Point& dir) + : + Function (f.n_components, f.get_time()), + f(f), + direction(dir) +{ + set_h(); + set_formula(); +} + + + +template +void +FunctionDerivative::set_h (double h_new) +{ + h = h_new; + incr = h*direction; +} + + + +template +void +FunctionDerivative::set_formula (DifferenceFormula form) +{ + formula = form; +} + + + +template +double +FunctionDerivative::value (const Point &p, + const unsigned int component) const +{ + switch (formula) + { + case Euler: + return (f.value(p+incr, component)-f.value(p-incr, component))/2*h; + case UpwindEuler: + return (f.value(p, component)-f.value(p-incr, component))/h; + default: + Assert(false, ExcInvalidFormula()); + } + return 0.; +} + + + +template +void +FunctionDerivative::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n = points.size(); + + switch (formula) + { + case Euler: + { + vector > p1(n); + vector > p2(n); + vector e2(n); + for (unsigned int i=0;i > p2(n); + vector e2(n); + for (unsigned int i=0;i; +template FunctionDerivative<2>; +template FunctionDerivative<3>; -- 2.39.5