From 9d7bf526e2fed05a73dd688c59c958e2043349d6 Mon Sep 17 00:00:00 2001 From: kanschat Date: Thu, 2 Aug 2007 15:11:31 +0000 Subject: [PATCH] Extended FunctionDerivative to compute gradients as well git-svn-id: https://svn.dealii.org/trunk@14881 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/base/auto_derivative_function.h | 25 ++- .../base/include/base/function_derivative.h | 38 ++-- deal.II/base/source/function_derivative.cc | 202 ++++++++++-------- deal.II/doc/news/changes.html | 17 ++ 4 files changed, 160 insertions(+), 122 deletions(-) diff --git a/deal.II/base/include/base/auto_derivative_function.h b/deal.II/base/include/base/auto_derivative_function.h index 66c8d93edf..e6fb5cc7d7 100644 --- a/deal.II/base/include/base/auto_derivative_function.h +++ b/deal.II/base/include/base/auto_derivative_function.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -82,15 +82,32 @@ class AutoDerivativeFunction : public Function enum DifferenceFormula { /** - * Symmetric Euler scheme + * The symmetric Euler + * formula of second order: + * @f[ + * u'(t) \approx + * \frac{u(t+h) - + * u(t-h)}{2h}. + * @f] */ Euler, /** - * Upwind Euler scheme + * The upwind Euler + * formula of first order: + * @f[ + * u'(t) \approx + * \frac{u(t) - + * u(t-h)}{h}. + * @f] */ UpwindEuler, /** - * Difference formula of 4th order. + * The fourth order scheme + * @f[ + * u'(t) \approx + * \frac{u(t-2h) - 8u(t-h) + * + 8u(t+h) - u(t+2h)}{12h}. + * @f] */ FourthOrder }; diff --git a/deal.II/base/include/base/function_derivative.h b/deal.II/base/include/base/function_derivative.h index feab8e9133..6a253ba425 100644 --- a/deal.II/base/include/base/function_derivative.h +++ b/deal.II/base/include/base/function_derivative.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -16,6 +16,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -37,20 +38,9 @@ DEAL_II_NAMESPACE_OPEN * @author Guido Kanschat, 2000 */ template -class FunctionDerivative : public Function +class FunctionDerivative : public AutoDerivativeFunction { public: - /** - * Names of difference formulas. - */ - enum DifferenceFormula - { - Euler, - UpwindEuler, - FourthOrder - }; - - /** * Constructor. Provided are the * functions to compute @@ -100,18 +90,20 @@ class FunctionDerivative : public Function * fourth order formula * (FourthOrder). */ - void set_formula (DifferenceFormula formula = Euler); - + void set_formula (typename AutoDerivativeFunction::DifferenceFormula formula + = AutoDerivativeFunction::Euler); /** - * Function value at one point. + * Change the base step size of + * the difference formula */ - virtual double value (const Point &p, - const unsigned int component = 0) const; + void set_h (const double h); + + virtual double value (const Point& p, + const unsigned int component = 0) const; + + virtual void vector_value(const Point& p, + Vector& value) const; - /** - * Function values at multiple - * points. - */ virtual void value_list (const std::vector > &points, std::vector &values, const unsigned int component = 0) const; @@ -157,7 +149,7 @@ class FunctionDerivative : public Function /** * Difference formula. */ - DifferenceFormula formula; + typename AutoDerivativeFunction::DifferenceFormula formula; /** * Helper object. Contains the diff --git a/deal.II/base/source/function_derivative.cc b/deal.II/base/source/function_derivative.cc index 3959e0c85f..0c9fe0ae84 100644 --- a/deal.II/base/source/function_derivative.cc +++ b/deal.II/base/source/function_derivative.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003, 2005, 2006 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -14,6 +14,7 @@ #include #include +#include #include @@ -24,7 +25,7 @@ FunctionDerivative::FunctionDerivative (const Function &f, const Point &dir, const double h) : - Function (f.n_components, f.get_time()), + AutoDerivativeFunction (h, f.n_components, f.get_time()), f(f), h(h), incr(1, h*dir) @@ -39,7 +40,7 @@ FunctionDerivative::FunctionDerivative (const Function& f, const std::vector >& dir, const double h) : - Function (f.n_components, f.get_time()), + AutoDerivativeFunction (h, f.n_components, f.get_time()), f(f), h(h), incr(dir.size()) @@ -53,28 +54,39 @@ FunctionDerivative::FunctionDerivative (const Function& f, template void -FunctionDerivative::set_formula (DifferenceFormula form) +FunctionDerivative::set_formula (typename AutoDerivativeFunction::DifferenceFormula form) { formula = form; } +template +void +FunctionDerivative::set_h (const double new_h) +{ + for (unsigned int i=0;i double FunctionDerivative::value (const Point &p, const unsigned int component) const { Assert (incr.size() == 1, - ExcMessage ("FunctionDerivative was not initialized for constant direection")); + ExcMessage ("FunctionDerivative was not initialized for constant direction")); switch (formula) { - case Euler: + case AutoDerivativeFunction::Euler: return (f.value(p+incr[0], component)-f.value(p-incr[0], component))/(2*h); - case UpwindEuler: + case AutoDerivativeFunction::UpwindEuler: return (f.value(p, component)-f.value(p-incr[0], component))/h; - case FourthOrder: + case AutoDerivativeFunction::FourthOrder: 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: @@ -85,6 +97,48 @@ FunctionDerivative::value (const Point &p, +template +void +FunctionDerivative::vector_value ( + const Point &p, + Vector& result) const +{ + Assert (incr.size() == 1, + ExcMessage ("FunctionDerivative was not initialized for constant direction")); + Vector aux(result.size()); + + // Formulas are the same as in + // value, but here we have to use + // Vector arithmetic + switch (formula) + { + case AutoDerivativeFunction::Euler: + f.vector_value(p+incr[0], result); + f.vector_value(p-incr[0], aux); + result.sadd(1./(2*h), -1./(2*h), aux); + return; + case AutoDerivativeFunction::UpwindEuler: + f.vector_value(p, result); + f.vector_value(p-incr[0], aux); + result.sadd(1./h, -1./h, aux); + return; + case AutoDerivativeFunction::FourthOrder: + f.vector_value(p-2*incr[0], result); + f.vector_value(p+2*incr[0], aux); + result.add(-1., aux); + f.vector_value(p-incr[0], aux); + result.add(-8., aux); + f.vector_value(p+incr[0], aux); + result.add(8., aux); + result.scale(1./(12*h)); + return; + default: + Assert(false, ExcInvalidFormula()); + } +} + + + template void FunctionDerivative::value_list (const std::vector > &points, @@ -97,95 +151,53 @@ FunctionDerivative::value_list (const std::vector > &points, Assert (incr.size() == points.size(), ExcDimensionMismatch(incr.size(), points.size())); + // Vector of auxiliary values + std::vector aux(n); + // Vector of auxiliary points + std::vector > paux(n); + // Use the same formulas as in + // value, but with vector + // arithmetic switch (formula) { - case Euler: - { - // let p1 and p2 be arrays of - // evaluation points shifted - // a little in direction j - std::vector > p1 = points; - std::vector > p2 = points; - - for (unsigned int i=0; i values2(n); - f.value_list(p1, values, component); - f.value_list(p2, values2, component); - - // finally compute finite - // differences - for (unsigned int i=0; i > p2 = points; - for (unsigned int i=0; i values2(n); - f.value_list(points, values, component); - f.value_list(p2, values2, component); - - // compute finite differences - for (unsigned int i=0; i > p_p = points; - std::vector > p_pp(n); - std::vector > p_m = points; - std::vector > p_mm(n); - for (unsigned int i=0;i e_p(n); - std::vector e_pp(n); - std::vector e_m(n); - - f.value_list(p_mm, values, component); - f.value_list(p_pp, e_pp, component); - f.value_list(p_p, e_p, component); - f.value_list(p_m, e_m, component); - - // compute finite differences - for (unsigned int i=0; i::Euler: + for (unsigned int i=0; i::UpwindEuler: + f.value_list(points, values, component); + for (unsigned int i=0; i::FourthOrder: + for (unsigned int i=0; i +

  • Changed: Implementing gradients for the class FunctionDerivative, it became evident that its enums for + difference formulas clashed with those of AutoDerivativeFunction. Therfore, only the latter + survived. +
    + (GK 2007/08/02) +

    +

  • Changed: When new multigrid transfer classes were introduced, the existing class MGTransferSelect was moved to the new header file @@ -437,6 +446,14 @@ inconvenience this causes.

    base

      + +

    1. Improved: FunctionDerivative is now + derived from AutoDerivativeFunction and implements + gradients as well, giving you automatic second derivatives of a function. +
      + (GK 2007/08/02) +

      +

    2. New: The function Utilities::fixed_power<n>(q) calculates q to the power of n where n is a constant known at compile time. It allows to -- 2.39.5