FunctionDerivative: probably buggy

git-svn-id: https://svn.dealii.org/trunk@3479 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
new file mode 100644 (file)
index 0000000..e080a5a
--- /dev/null
+++ b/deal.II/base/include/base/function_derivative.h
@@ -0,0 +1,111 @@
+//-------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------------
+#ifndef __deal2__function_derivative_h
+#define __deal2__function_derivative_h
+
+#include <base/exceptions.h>
+#include <base/function.h>
+
+/**
+ * Names of difference formulas.
+ */
+enum DifferenceFormula
+{
+  Euler,
+  UpwindEuler
+};
+
+
+/**
+ * Derivative of a function object.  The value access functions of
+ * this class return the directional derivative of a function with
+ * respect to a direction provided on construction. If @p{b} is the
+ * vector, the derivative @p{b . grad f} is computed. This derivative
+ * is evaluated directly, not by computing the gradient of @p{f} and
+ * its scalar product with @p{b}.
+ *
+ * The derivative is computed numerically, using one of the provided
+ * difference formulas.
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <int dim>
+class FunctionDerivative : public Function<dim>
+{
+public:
+                                  /**
+                                   * Constructor. Provided are the
+                                   * function to compute derivatives
+                                   * of and the direction vector of
+                                   * the differentiation.
+                                   */
+  FunctionDerivative (const Function<dim>& f,
+                     const Point<dim>& direction);
+
+                                  /**
+                                   * Set step size of the difference
+                                   * formula. This is set to the
+                                   * default in the constructor.
+                                   */
+  void set_h (double h_new = 1.e-6);
+  
+                                  /**
+                                   * Choose the difference formula.
+                                   * This is set to the default in
+                                   * the constructor.
+                                   */
+  void set_formula (DifferenceFormula formula = Euler);
+  
+                                  /**
+                                   * Function value at one point.
+                                   */
+  virtual double value (const Point<dim>   &p,
+                       const unsigned int  component = 0) const;
+  
+                                  /**
+                                   * Function values at multiple points.
+                                   */
+  virtual void value_list (const vector<Point<dim> > &points,
+                          vector<double>            &values,
+                          const unsigned int         component = 0) const;
+
+
+  DeclException0(ExcInvalidFormula);
+  
+private:
+                                  /**
+                                   * Function for differentiation.
+                                   */
+  const Function<dim>& f;
+                                  /**
+                                   * Differentiation vector.
+                                   */
+  Point<dim> direction;
+  
+                                  /**
+                                   * Step size of the difference formula.
+                                   */
+  double h;
+                                  /**
+                                   * Difference formula.
+                                   */
+  DifferenceFormula formula;
+                                  /**
+                                   * Helper object. Contains the
+                                   * increment vector for the
+                                   * formula.
+                                   */
+  Point<dim> incr;
+};
+
+#endif

diff --git a/deal.II/base/source/function_derivative.cc b/deal.II/base/source/function_derivative.cc
new file mode 100644 (file)
index 0000000..5866eb9
--- /dev/null
+++ b/deal.II/base/source/function_derivative.cc
@@ -0,0 +1,118 @@
+//--------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//--------------------------------------------------------------------------
+
+
+#include <base/point.h>
+#include <base/function_derivative.h>
+
+#include <cmath>
+
+template <int dim>
+FunctionDerivative<dim>::FunctionDerivative (const Function<dim>& f,
+                                            const Point<dim>& dir)
+  :
+  Function<dim> (f.n_components, f.get_time()),
+  f(f),
+  direction(dir)
+{
+  set_h();
+  set_formula();
+}
+
+
+
+template <int dim>
+void
+FunctionDerivative<dim>::set_h (double h_new)
+{
+  h = h_new;
+  incr = h*direction;
+}
+
+
+
+template <int dim>
+void
+FunctionDerivative<dim>::set_formula (DifferenceFormula form)
+{
+  formula = form;
+}
+
+
+
+template <int dim>
+double
+FunctionDerivative<dim>::value (const Point<dim>   &p,
+                               const unsigned int  component) const
+{
+  switch (formula)
+    {
+    case Euler:
+      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;
+    default:
+      Assert(false, ExcInvalidFormula());
+    }
+  return 0.;
+}
+
+
+
+template <int dim>
+void
+FunctionDerivative<dim>::value_list (const vector<Point<dim> > &points,
+                                    vector<double>            &values,
+                                    const unsigned int         component) const
+{
+  const unsigned int n = points.size();
+  
+  switch (formula)
+    {
+    case Euler:
+    {
+      vector<Point<dim> > p1(n);
+      vector<Point<dim> > p2(n);
+      vector<double> e2(n);
+      for (unsigned int i=0;i<n;++i)
+       {
+         p1[i] = points[i]+incr;
+         p2[i] = points[i]-incr;
+       }
+      f.value_list(p1, values, component);
+      f.value_list(p2, e2, component);
+      for (unsigned int i=0;i<n;++i)
+       values[i] = (values[i]-e2[i])/2*h;
+      break;
+    }    
+    case UpwindEuler:
+    {
+      vector<Point<dim> > p2(n);
+      vector<double> e2(n);
+      for (unsigned int i=0;i<n;++i)
+       p2[i] = points[i]-incr;
+      f.value_list(points, values, component);
+      f.value_list(p2, e2, component);
+      for (unsigned int i=0;i<n;++i)
+       values[i] = (values[i]-e2[i])/h;
+      break;
+    }
+    default:
+      Assert(false, ExcInvalidFormula());
+    }
+}
+
+
+template FunctionDerivative<1>;
+template FunctionDerivative<2>;
+template FunctionDerivative<3>;