From: Guido Kanschat Date: Mon, 13 Nov 2000 14:32:44 +0000 (+0000) Subject: 4th order formula in FunctionDerivative X-Git-Tag: v8.0.0~19947 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=faad29060b61f05345fda6bb1c03481c3e195d0b;p=dealii.git 4th order formula in FunctionDerivative git-svn-id: https://svn.dealii.org/trunk@3485 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/function_derivative.h b/deal.II/base/include/base/function_derivative.h index e080a5ac4c..2d0f7aca93 100644 --- a/deal.II/base/include/base/function_derivative.h +++ b/deal.II/base/include/base/function_derivative.h @@ -22,7 +22,8 @@ enum DifferenceFormula { Euler, - UpwindEuler + UpwindEuler, + FourthOrder }; @@ -35,7 +36,9 @@ enum DifferenceFormula * its scalar product with @p{b}. * * The derivative is computed numerically, using one of the provided - * difference formulas. + * difference formulas (see @p{set_formula} for available + * schemes). Experimenting with @p{h} and the difference scheme may be + * necessary to obtain sufficient results. * * @author Guido Kanschat, 2000 */ @@ -63,6 +66,13 @@ public: * Choose the difference formula. * This is set to the default in * the constructor. + * + * Formulas implemented right now + * are first order backward Euler + * (@p{UpwindEuler}), second order + * symmetric Euler (@p{Euler}) and + * a symmetric fourth order formula + * (@p{FourthOrder}). */ void set_formula (DifferenceFormula formula = Euler); diff --git a/deal.II/base/source/function_derivative.cc b/deal.II/base/source/function_derivative.cc index ad76305d02..5a4636de58 100644 --- a/deal.II/base/source/function_derivative.cc +++ b/deal.II/base/source/function_derivative.cc @@ -48,7 +48,7 @@ FunctionDerivative::set_formula (DifferenceFormula form) formula = form; } - +//TODO: Discussion on an efficient implementation of Point additions. template double @@ -61,6 +61,9 @@ FunctionDerivative::value (const Point &p, 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; + 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); default: Assert(false, ExcInvalidFormula()); } @@ -110,6 +113,34 @@ FunctionDerivative::value_list (const vector > &points, values[i] = (values[i]-e2[i])/h; break; } + case FourthOrder: + { + vector > p_p(n); + vector > p_pp(n); + vector > p_m(n); + vector > p_mm(n); + vector e_p(n); + vector e_pp(n); + vector e_m(n); + for (unsigned int i=0;i