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 e080a5ac4c907d9af3a0012bb83506b365a23aaf..2d0f7aca93d1d5b61a6ca756740954e77a9d3eac 100644 (file)
--- 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 ad76305d02fd0eacf646257afc238660e50a3605..5a4636de585ae641c67eafe318bed3767d2368ff 100644 (file)
--- 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());
     }