From 269fc530f4ba13cf7f5c646531b21e5b9ea2c287 Mon Sep 17 00:00:00 2001
From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 13 Nov 2000 14:32:44 +0000
Subject: [PATCH] 4th order formula in FunctionDerivative

git-svn-id: https://svn.dealii.org/trunk@3485 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../base/include/base/function_derivative.h   | 14 ++++++--
 deal.II/base/source/function_derivative.cc    | 33 ++++++++++++++++++-
 2 files changed, 44 insertions(+), 3 deletions(-)

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<dim>::set_formula (DifferenceFormula form)
   formula = form;
 }
 
-
+//TODO: Discussion on an efficient implementation of Point additions.
 
 template <int dim>
 double
@@ -61,6 +61,9 @@ FunctionDerivative<dim>::value (const Point<dim>   &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<dim>::value_list (const vector<Point<dim> > &points,
 	values[i] = (values[i]-e2[i])/h;
       break;
     }
+    case FourthOrder:
+    {
+      vector<Point<dim> > p_p(n);
+      vector<Point<dim> > p_pp(n);
+      vector<Point<dim> > p_m(n);
+      vector<Point<dim> > p_mm(n);
+      vector<double> e_p(n);
+      vector<double> e_pp(n);
+      vector<double> e_m(n);
+      for (unsigned int i=0;i<n;++i)
+	{
+	  p_p[i] = points[i]+incr;
+	  p_pp[i] = p_p[i]+incr;
+	  p_m[i] = points[i]-incr;
+	  p_mm[i] = p_m[i]-incr;
+	}
+      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);
+      
+      for (unsigned int i=0;i<n;++i)
+	{
+	  values[i] = (values[i]-p_pp[i]+8*(p_p[i]-p_m[i]))/(12*h);
+	}
+      break;
+    }    
+
     default:
       Assert(false, ExcInvalidFormula());
     }
-- 
2.39.5