From 820de33707c0e75325e519d9d9f4ff22426df70f Mon Sep 17 00:00:00 2001
From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 12 Nov 2000 21:36:25 +0000
Subject: [PATCH] FunctionDerivative: probably buggy

git-svn-id: https://svn.dealii.org/trunk@3479 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../base/include/base/function_derivative.h   | 111 ++++++++++++++++
 deal.II/base/source/function_derivative.cc    | 118 ++++++++++++++++++
 2 files changed, 229 insertions(+)
 create mode 100644 deal.II/base/include/base/function_derivative.h
 create mode 100644 deal.II/base/source/function_derivative.cc

diff --git a/deal.II/base/include/base/function_derivative.h b/deal.II/base/include/base/function_derivative.h
new file mode 100644
index 0000000000..e080a5ac4c
--- /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
index 0000000000..5866eb9747
--- /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>;
-- 
2.39.5