Extended FunctionDerivative to compute gradients as well

git-svn-id: https://svn.dealii.org/trunk@14881 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/base/include/base/auto_derivative_function.h b/deal.II/base/include/base/auto_derivative_function.h
index 66c8d93edfedd3687930d5b2be4594957452a057..e6fb5cc7d7200f9c17ded1c8b598dcf3da363f20 100644 (file)
--- 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<dim>
     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 feab8e9133646821ea7866a20838caa6255b8c7b..6a253ba4256214d0b075edbc6213194554d2defd 100644 (file)
--- 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 <base/config.h>
 #include <base/exceptions.h>
 #include <base/function.h>
+#include <base/auto_derivative_function.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -37,20 +38,9 @@ DEAL_II_NAMESPACE_OPEN
  * @author Guido Kanschat, 2000
  */
 template <int dim>
-class FunctionDerivative : public Function<dim>
+class FunctionDerivative : public AutoDerivativeFunction<dim>
 {
   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<dim>
                                      * fourth order formula
                                      * (<tt>FourthOrder</tt>).
                                      */
-    void set_formula (DifferenceFormula formula = Euler);
-    
+    void set_formula (typename AutoDerivativeFunction<dim>::DifferenceFormula formula
+                     = AutoDerivativeFunction<dim>::Euler);
                                     /**
-                                     * Function value at one point.
+                                     * Change the base step size of
+                                     * the difference formula
                                      */
-    virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
+    void set_h (const double h);
+    
+    virtual double value (const Point<dim>& p,
+                         const unsigned int component = 0) const;
+    
+    virtual void vector_value(const Point<dim>& p,
+                             Vector<double>& value) const;
     
-                                    /**
-                                     * Function values at multiple
-                                     * points.
-                                     */
     virtual void value_list (const std::vector<Point<dim> > &points,
                             std::vector<double>            &values,
                             const unsigned int              component = 0) const;    
@@ -157,7 +149,7 @@ class FunctionDerivative : public Function<dim>
                                     /**
                                      * Difference formula.
                                      */
-    DifferenceFormula formula;
+    typename AutoDerivativeFunction<dim>::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 3959e0c85fee693801b0a82c727c480f840ce093..0c9fe0ae8455dade7bee56403fefaddc4ce7d172 100644 (file)
--- 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 <base/point.h>
 #include <base/function_derivative.h>
+#include <lac/vector.h>
 
 #include <cmath>
 
@@ -24,7 +25,7 @@ FunctionDerivative<dim>::FunctionDerivative (const Function<dim> &f,
                                             const Point<dim>    &dir,
                                             const double         h)
                :
-               Function<dim> (f.n_components, f.get_time()),
+               AutoDerivativeFunction<dim> (h, f.n_components, f.get_time()),
                 f(f),
                 h(h),
                 incr(1, h*dir)
@@ -39,7 +40,7 @@ FunctionDerivative<dim>::FunctionDerivative (const Function<dim>& f,
                                             const std::vector<Point<dim> >& dir,
                                             const double h)
                :
-               Function<dim> (f.n_components, f.get_time()),
+               AutoDerivativeFunction<dim> (h, f.n_components, f.get_time()),
   f(f),
   h(h),
   incr(dir.size())
@@ -53,28 +54,39 @@ FunctionDerivative<dim>::FunctionDerivative (const Function<dim>& f,
 
 template <int dim>
 void
-FunctionDerivative<dim>::set_formula (DifferenceFormula form)
+FunctionDerivative<dim>::set_formula (typename AutoDerivativeFunction<dim>::DifferenceFormula form)
 {
   formula = form;
 }
 
 
 
+template <int dim>
+void
+FunctionDerivative<dim>::set_h (const double new_h)
+{
+  for (unsigned int i=0;i<incr.size ();++i)
+    incr[i] *= new_h/h;
+  h = new_h;
+}
+
+
+
 template <int dim>
 double
 FunctionDerivative<dim>::value (const Point<dim>   &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<dim>::Euler:
            return (f.value(p+incr[0], component)-f.value(p-incr[0], component))/(2*h);
-      case UpwindEuler:
+      case AutoDerivativeFunction<dim>::UpwindEuler:
            return (f.value(p, component)-f.value(p-incr[0], component))/h;
-      case FourthOrder:
+      case AutoDerivativeFunction<dim>::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<dim>::value (const Point<dim>   &p,
 
 
 
+template <int dim>
+void
+FunctionDerivative<dim>::vector_value (
+  const Point<dim>   &p,
+  Vector<double>& result) const
+{
+  Assert (incr.size() == 1,
+         ExcMessage ("FunctionDerivative was not initialized for constant direction"));
+  Vector<double> aux(result.size());
+
+                                  // Formulas are the same as in
+                                  // value, but here we have to use
+                                  // Vector arithmetic
+  switch (formula)
+    {
+      case AutoDerivativeFunction<dim>::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<dim>::UpwindEuler:
+           f.vector_value(p, result);
+           f.vector_value(p-incr[0], aux);
+           result.sadd(1./h, -1./h, aux);
+           return;
+      case AutoDerivativeFunction<dim>::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 <int dim>
 void
 FunctionDerivative<dim>::value_list (const std::vector<Point<dim> > &points,
@@ -97,95 +151,53 @@ FunctionDerivative<dim>::value_list (const std::vector<Point<dim> > &points,
     Assert (incr.size() == points.size(),
            ExcDimensionMismatch(incr.size(), points.size()));
 
+                                  // Vector of auxiliary values
+  std::vector<double> aux(n);
+                                  // Vector of auxiliary points
+  std::vector<Point<dim> > 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<Point<dim> > p1 = points;
-       std::vector<Point<dim> > p2 = points;
-
-       for (unsigned int i=0; i<n; ++i)
-         {
-           const unsigned int j = (variable_direction) ? i:0;
-           p1[i] += incr[j];
-           p2[i] -= incr[j];
-         }     
-
-                                        // next get values of
-                                        // functions at these points
-       std::vector<double> values2(n);
-       f.value_list(p1, values,  component);
-       f.value_list(p2, values2, component);
-
-                                        // finally compute finite
-                                        // differences
-       for (unsigned int i=0; i<n; ++i)
-         values[i] = (values[i]-values2[i])/(2*h);
-         
-       break;
-      }
-       
-      case UpwindEuler:
-      {
-                                        // compute upwind points
-       std::vector<Point<dim> > p2 = points;
-       for (unsigned int i=0; i<n; ++i)
-         {
-           const unsigned int j = (variable_direction) ? i:0;
-           p2[i] -= incr[j];
-         }
-
-                                        // get values at points
-       std::vector<double> values2(n);
-       f.value_list(points, values,  component);
-       f.value_list(p2,     values2, component);
-
-                                        // compute finite differences
-       for (unsigned int i=0; i<n; ++i)
-         values[i] = (values[i]-values2[i])/h;
-       break;
-      }
-       
-      case FourthOrder:
-      {
-                                        // first compute evaluation
-                                        // points
-       std::vector<Point<dim> > p_p = points;
-       std::vector<Point<dim> > p_pp(n);
-       std::vector<Point<dim> > p_m = points;
-       std::vector<Point<dim> > p_mm(n);
-       for (unsigned int i=0;i<n;++i)
-         {
-           const unsigned int j = (variable_direction) ? i:0;
-           p_p[i] += incr[j];
-           p_pp[i] = p_p[i]+incr[j];
-           p_m[i] -= incr[j];
-           p_mm[i] = p_m[i]-incr[j];
-         }
-
-                                        // next compute values of
-                                        // function at these
-                                        // points. use @p{values} for
-                                        // @p{e_mm}
-       std::vector<double> e_p(n);
-       std::vector<double> e_pp(n);
-       std::vector<double> 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<n; ++i)
-         values[i] = (values[i]-e_pp[i]+8*(e_p[i]-e_m[i]))/(12*h);
-         
-       break;
-      }    
-
+      case AutoDerivativeFunction<dim>::Euler:
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]+incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, values, component);
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]-incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, aux, component);
+           for (unsigned int i=0; i<n; ++i)
+             values[i] = (values[i]-aux[i])/(2*h);
+           return;
+      case AutoDerivativeFunction<dim>::UpwindEuler:
+           f.value_list(points, values, component);
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]-incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, aux, component);
+           for (unsigned int i=0; i<n; ++i)
+             values[i] = (values[i]-aux[i])/h;
+           return;
+      case AutoDerivativeFunction<dim>::FourthOrder:
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]-2*incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, values, component);
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]+2*incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, aux, component);
+           for (unsigned int i=0; i<n; ++i)
+             values[i] -= aux[i];
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]+incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, aux, component);
+           for (unsigned int i=0; i<n; ++i)
+             values[i] += 8.*aux[i];
+           for (unsigned int i=0; i<n; ++i)
+             paux[i] = points[i]-incr[(variable_direction) ? i : 0U];
+           f.value_list(paux, aux, component);
+           for (unsigned int i=0; i<n; ++i)
+             values[i] = (values[i] - 8.*aux[i])/(12*h);
+           return;
       default:
            Assert(false, ExcInvalidFormula());
     }

diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html
index 48d5cf736eae9b55747309fab9e4fe3053fc7d28..57a95adebd437d220719fdd46ee5af2a7cc7f8ce 100644 (file)
--- a/deal.II/doc/news/changes.html
+++ b/deal.II/doc/news/changes.html
@@ -45,6 +45,15 @@ inconvenience this causes.
 
 <ol>
 
+  <li> <p>Changed: Implementing gradients for the class <code
+  class="class">FunctionDerivative</code>, it became evident that its enums for
+  difference formulas clashed with those of <code
+  class="class">AutoDerivativeFunction</code>. Therfore, only the latter
+  survived.
+  <br>
+  (GK 2007/08/02)
+  </p>
+
   <li> <p>Changed: When new multigrid transfer classes were introduced,
   the existing class <code class="class">MGTransferSelect</code> was
   moved to the new header file
@@ -437,6 +446,14 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+
+  <li> <p> Improved: <code class="class">FunctionDerivative</code> is now
+  derived from <code class="class">AutoDerivativeFunction</code> and implements
+  gradients as well, giving you automatic second derivatives of a function.
+  <br>
+  (GK 2007/08/02)
+  </p>
+
   <li> <p> New: The function <code>Utilities::fixed_power&lt;n&gt;(q)</code>
        calculates <code>q</code> to the power of <code>n</code> where
        <code>n</code> is a constant known at compile time. It allows to