]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implementation of new AutoDerivativeFunction.
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Tue, 15 May 2001 16:08:21 +0000 (16:08 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Tue, 15 May 2001 16:08:21 +0000 (16:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@4610 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/auto_derivative_function.h [new file with mode: 0644]
deal.II/base/source/auto_derivative_function.cc [new file with mode: 0644]

diff --git a/deal.II/base/include/base/auto_derivative_function.h b/deal.II/base/include/base/auto_derivative_function.h
new file mode 100644 (file)
index 0000000..91137b9
--- /dev/null
@@ -0,0 +1,257 @@
+//----------------------------  auto_derivative_function.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal 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.
+//
+//----------------------------  auto_derivative_function.h  ---------------------------
+#ifndef __deal2__auto_derivative_function_h
+#define __deal2__auto_derivative_function_h
+
+
+#include <base/exceptions.h>
+#include <base/function.h>
+
+
+/**
+ * This class automatically computes the gradient of a function by
+ * employing numerical difference quotients. This only, if the user
+ * function does not provide the gradient function himself.
+ *
+ * @sect3{Usage}
+ * The following example of an user defined function overloads and
+ * implements only the @p{value} function but not the @p{gradient}
+ * function. If the @p{gradient} function is invoked then the gradient
+ * function implemented by the @p{AutoDerivativeFunction} is called,
+ * where the latter function imployes numerical difference quotients.
+ *
+ * @begin{verbatim}
+ * class UserFunction: public AutoDerivativeFunction
+ * {               // access to one component at one point
+ *   double value (const Point<dim> &p, const
+ *                 unsigned int component = 0) const
+ *          { // Implementation ....  };
+ * } user_function;
+ *
+ *            // gradient by employing difference quotients.
+ * Tensor<1,dim> grad=user_function.gradient(some_point);
+ * @end{verbatim}
+ * 
+ * If the user overloads and implements also the gradient function,
+ * then, of course, the users gradient function is called.
+ *
+ * Note, that the usage of the @p{value} and @p{gradient} functions
+ * explained above, also applies to the @p{value_list} and
+ * @p{gradient_list} functions as well as to the vector valued
+ * versions of these functions, see e.g. @p{vector_value},
+ * @p{vector_gradient}, @p{vector_value_list} and
+ * @p{vector_gradient_list}.
+ *
+ * The @p{gradient} and @p{gradient_list} functions make use of the
+ * @p{value} function. The @p{vector_gradient} and
+ * @p{vector_gradient_list} make use of the @p{vector_value}
+ * function. Make sure that the user defined function implements the
+ * @p{value} function and the @p{vector_value} function, respectively.
+ *
+ * Furthermore note, that an object of this class does @em{not} represent
+ * the derivative of a function, like @ref{FunctionDerivative}, that
+ * gives a directional derivate by calling the @p{value} function. In
+ * fact, this class (the @p{AutoDerivativeFunction} class) can
+ * substitute the @p{Function} class as base class for user defined
+ * classes. This class implements the @p{gradient} functions for
+ * automatic computation of numerical difference quotients and serves
+ * as intermediate class between the base @p{Function} class and the
+ * user defined function class.
+ *
+ * @author Ralf Hartmann, 2001
+ */
+template <int dim>
+class AutoDerivativeFunction : public Function<dim>
+{
+  public:
+    
+                                    /**
+                                     * Names of difference formulas.
+                                     */
+    enum DifferenceFormula
+    {
+         Euler,
+         UpwindEuler,
+         FourthOrder
+    };
+
+                                    /**
+                                     * Constructor. Takes the
+                                     * difference step size
+                                     * @p{h}. It's within the user's
+                                     * responsibility to choose an
+                                     * appropriate value here. @p{h}
+                                     * should be chosen taking into
+                                     * account the absolute value as
+                                     * well as the amount of local
+                                     * variation of the function.
+                                     * Setting @p{h=1e-6} might be a
+                                     * good choice for functions with
+                                     * an absolute value of about 1,
+                                     * that furthermore does not vary
+                                     * to much.
+                                     *
+                                     * @p{h} can be changed later
+                                     * using the @p{set_h} function.
+                                     *
+                                     * Sets @p{DifferenceFormula}
+                                     * @p{formula} to the default
+                                     * @p{Euler} formula of the
+                                     * @p{set_formula}
+                                     * function. Change this preset
+                                     * formula by calling the
+                                     * @p{set_formula} function.
+                                     */
+    AutoDerivativeFunction (const double h,
+                           const unsigned int n_components = 1,
+                           const double       initial_time = 0.0);
+    
+                                    /**
+                                     * Virtual destructor; absolutely
+                                     * necessary in this case.
+                                     */
+    virtual ~AutoDerivativeFunction ();
+    
+                                    /**
+                                     * Choose the difference formula.
+                                     *
+                                     * 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 (const DifferenceFormula formula = Euler);
+
+                                    /**
+                                     * Takes the difference step size
+                                     * @p{h}. It's within the user's
+                                     * responsibility to choose an
+                                     * appropriate value here. @p{h}
+                                     * should be chosen taking into
+                                     * account the absolute value of
+                                     * as well as the amount of local
+                                     * variation of the function.
+                                     * Setting @p{h=1e-6} might be a
+                                     * good choice for functions with
+                                     * an absolute value of about 1,
+                                     * that furthermore does not vary
+                                     * to much.
+                                     */
+    void set_h (const double h);
+    
+                                    /**
+                                     * Return the gradient of the
+                                     * specified component of the
+                                     * function at the given point.
+                                     *
+                                     * Imployes numerical difference
+                                     * quotients using the preset
+                                     * @p{DifferenceFormula}
+                                     * @p{formula}.
+                                     */
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                   const unsigned int  component = 0) const;
+
+                                    /**
+                                     * Return the gradient of all
+                                     * components of the
+                                     * function at the given point.
+                                     *
+                                     * Imployes numerical difference
+                                     * quotients using the preset
+                                     * @p{DifferenceFormula}
+                                     * @p{formula}.
+                                     */
+    virtual void vector_gradient (const Point<dim>            &p,
+                                 typename std::vector<Tensor<1,dim> > &gradients) const;
+    
+                                    /**
+                                     * Set @p{gradients} to the
+                                     * gradients of the specified
+                                     * component of the function at
+                                     * the @p{points}.  It is assumed
+                                     * that @p{gradients} already has the
+                                     * right size, i.e.  the same
+                                     * size as the @p{points} array.
+                                     *
+                                     * Imployes numerical difference
+                                     * quotients using the preset
+                                     * @p{DifferenceFormula}
+                                     * @p{formula}.
+                                     */
+    virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+                               typename std::vector<Tensor<1,dim> >    &gradients,
+                               const unsigned int              component = 0) const;
+    
+                                    /**
+                                     * Set @p{gradients} to the gradients of
+                                     * the function at the @p{points},
+                                     * for all components.
+                                     * It is assumed that @p{gradients} 
+                                     * already has the right size, i.e.
+                                     * the same size as the @p{points} array.
+                                     *
+                                     * The outer loop over
+                                     * @p{gradients} is over the points
+                                     * in the list, the inner loop
+                                     * over the different components
+                                     * of the function.
+                                     *
+                                     * Imployes numerical difference
+                                     * quotients using the preset
+                                     * @p{DifferenceFormula}
+                                     * @p{formula}.
+                                     */
+    virtual void vector_gradient_list (const typename std::vector<Point<dim> > &points,
+                                      typename std::vector<typename std::vector<Tensor<1,dim> > > &gradients) const;
+
+                                    /**
+                                     * Returns a
+                                     * @p{DifferenceFormula} of the
+                                     * order @p{ord} at minimum.
+                                     */
+    static DifferenceFormula get_formula_of_order(const unsigned int ord);
+
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException0(ExcInvalidFormula);
+
+  private:
+    
+                                    /**
+                                     * Step size of the difference
+                                     * formula. Set by the @p{set_h}
+                                     * function.
+                                     */
+    double h;
+
+                                    /**
+                                     * Includes the unit vectors
+                                     * scaled by @p{h}.
+                                     */
+    typename std::vector<Tensor<1,dim> > ht;
+    
+                                    /**
+                                     * Difference formula. Set by the
+                                     * @p{set_formula} function.
+                                     */
+    DifferenceFormula formula;
+};
+
+
+
+#endif
diff --git a/deal.II/base/source/auto_derivative_function.cc b/deal.II/base/source/auto_derivative_function.cc
new file mode 100644 (file)
index 0000000..59de4c2
--- /dev/null
@@ -0,0 +1,325 @@
+//--------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 by the deal 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/auto_derivative_function.h>
+#include <lac/vector.h>
+
+#include <cmath>
+
+
+// if necessary try to work around a bug in the IBM xlC compiler
+#ifdef XLC_WORK_AROUND_STD_BUG
+using namespace std;
+#endif
+
+
+template <int dim>
+AutoDerivativeFunction<dim>::AutoDerivativeFunction (const double hh,
+                                                    const unsigned int n_components,
+                                                    const double       initial_time):
+               Function<dim>(n_components, initial_time),
+               h(1),
+               ht(dim),
+               formula(Euler)
+{
+  set_h(hh);
+  set_formula();
+}
+
+
+template <int dim>
+AutoDerivativeFunction<dim>::~AutoDerivativeFunction ()
+{}
+
+
+
+template <int dim>
+void
+AutoDerivativeFunction<dim>::set_formula (const DifferenceFormula form)
+{
+  formula = form;
+}
+
+
+template <int dim>
+void
+AutoDerivativeFunction<dim>::set_h (const double hh)
+{
+  h=hh;
+  for (unsigned int i=0; i<dim; ++i)
+    ht[i][i]=h;
+}
+
+
+template <int dim>
+Tensor<1,dim>
+AutoDerivativeFunction<dim>::gradient (const Point<dim>   &p,
+                                      const unsigned int  comp) const
+{
+  Tensor<1,dim> grad;
+  switch (formula)
+    {
+      case UpwindEuler:
+      {
+       Tensor<1,dim> q1;
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           q1=p-ht[i];
+           grad[i]=(value(p, comp)-value(q1, comp))/h;
+         }
+       break;
+      }
+      case Euler:
+      {
+       Tensor<1,dim> q1, q2;
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           q1=p+ht[i];
+           q2=p-ht[i];
+           grad[i]=(value(q1, comp)-value(q2, comp))/(2*h);
+         }
+       break;
+      }
+      case FourthOrder:
+      {
+       Tensor<1,dim> q1, q2, q3, q4;
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           q2=p+ht[i];
+           q1=q2+ht[i];
+           q3=p-ht[i];
+           q4=q3-ht[i];
+           grad[i]=(-  value(q1, comp)
+                    +8*value(q2, comp)
+                    -8*value(q3, comp)
+                    +  value(q4, comp))/(12*h);
+           
+         }
+       break;
+      }
+      default:
+           Assert(false, ExcInvalidFormula());
+    }
+  return grad;
+}
+
+
+template <int dim>
+void AutoDerivativeFunction<dim>::vector_gradient (const Point<dim>       &p,
+                                                  typename std::vector<Tensor<1,dim> > &gradients) const
+{
+  Assert (gradients.size() == n_components, ExcDimensionMismatch(gradients.size(), n_components));
+  
+  switch (formula)
+    {
+      case UpwindEuler:
+      {
+       Tensor<1,dim> q1;
+       Vector<double> v(n_components), v1(n_components);
+       const double h_inv=1./h;
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           q1=p-ht[i];
+           vector_value(p, v);
+           vector_value(q1, v1);
+           
+           for (unsigned int comp=0; comp<n_components; ++comp)
+             gradients[comp][i]=(v(comp)-v1(comp))*h_inv;
+         }
+       break;
+      }
+      case Euler:
+      {
+       Tensor<1,dim> q1, q2;
+       Vector<double> v1(n_components), v2(n_components);
+       const double h_inv_2=1./(2*h);
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           q1=p+ht[i];
+           q2=p-ht[i];
+           vector_value(q1, v1);
+           vector_value(q2, v2);
+           
+           for (unsigned int comp=0; comp<n_components; ++comp)
+             gradients[comp][i]=(v1(comp)-v2(comp))*h_inv_2;
+         }
+       break;
+      }
+      case FourthOrder:
+      {
+       Tensor<1,dim> q1, q2, q3, q4;
+       Vector<double> v1(n_components), v2(n_components), v3(n_components), v4(n_components);
+       const double h_inv_12=1./(12*h);
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           q2=p+ht[i];
+           q1=q2+ht[i];
+           q3=p-ht[i];
+           q4=q3-ht[i];
+           vector_value(q1, v1);
+           vector_value(q2, v2);
+           vector_value(q3, v3);
+           vector_value(q4, v4);
+           
+           for (unsigned int comp=0; comp<n_components; ++comp)
+             gradients[comp][i]=(-v1(comp)+8*v2(comp)-8*v3(comp)+v4(comp))*h_inv_12;
+         }
+       break;
+      }
+      default:
+           Assert(false, ExcInvalidFormula());
+    }
+}
+
+
+template <int dim>
+void AutoDerivativeFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
+                                                typename std::vector<Tensor<1,dim> >    &gradients,
+                                                const unsigned int              comp) const
+{
+  Assert (gradients.size() == points.size(),
+         ExcDimensionMismatch(gradients.size(), points.size()));
+
+  switch (formula)
+    {
+      case UpwindEuler:
+      {
+       Tensor<1,dim> q1;
+       for (unsigned p=0; p<points.size(); ++p)
+         for (unsigned int i=0; i<dim; ++i)
+           {
+             q1=points[p]-ht[i];
+             gradients[p][i]=(value(points[p], comp)-value(q1, comp))/h;
+           }
+       break;
+      }
+      case Euler:
+      {
+       Tensor<1,dim> q1, q2;
+       for (unsigned p=0; p<points.size(); ++p)
+         for (unsigned int i=0; i<dim; ++i)
+           {
+             q1=points[p]+ht[i];
+             q2=points[p]-ht[i];
+             gradients[p][i]=(value(q1, comp)-value(q2, comp))/(2*h);
+           }
+       break;
+      }
+      case FourthOrder:
+      {
+       Tensor<1,dim> q1, q2, q3, q4;
+       for (unsigned p=0; p<points.size(); ++p)
+         for (unsigned int i=0; i<dim; ++i)
+           {
+             q2=points[p]+ht[i];
+             q1=q2+ht[i];
+             q3=points[p]-ht[i];
+             q4=q3-ht[i];
+             gradients[p][i]=(-  value(q1, comp)
+                              +8*value(q2, comp)
+                              -8*value(q3, comp)
+                              +  value(q4, comp))/(12*h);
+           }
+       break;
+      }
+      default:
+           Assert(false, ExcInvalidFormula());
+    }
+}
+
+
+
+template <int dim>
+void AutoDerivativeFunction<dim>::vector_gradient_list (const typename std::vector<Point<dim> >            &points,
+                                                       typename std::vector<std::vector<Tensor<1,dim> > > &gradients) const
+{
+  Assert (gradients.size() == points.size(),
+         ExcDimensionMismatch(gradients.size(), points.size()));
+  for (unsigned p=0; p<points.size(); ++p)
+    Assert (gradients[p].size() == n_components,
+           ExcDimensionMismatch(gradients.size(), n_components));
+    
+  switch (formula)
+    {
+      case UpwindEuler:
+      {
+       Tensor<1,dim> q1;
+       for (unsigned p=0; p<points.size(); ++p)
+         for (unsigned int i=0; i<dim; ++i)
+           {
+             q1=points[p]-ht[i];
+             for (unsigned int comp=0; comp<n_components; ++comp)
+               gradients[p][comp][i]=(value(points[p], comp)-value(q1, comp))/h;
+         }
+       break;
+      }
+      case Euler:
+      {
+       Tensor<1,dim> q1, q2;
+       for (unsigned p=0; p<points.size(); ++p)
+         for (unsigned int i=0; i<dim; ++i)
+           {
+             q1=points[p]+ht[i];
+             q2=points[p]-ht[i];
+             for (unsigned int comp=0; comp<n_components; ++comp)
+               gradients[p][comp][i]=(value(q1, comp)-value(q2, comp))/(2*h);
+           }
+       break;
+      }
+      case FourthOrder:
+      {
+       Tensor<1,dim> q1, q2, q3, q4;
+       for (unsigned p=0; p<points.size(); ++p)
+         for (unsigned int i=0; i<dim; ++i)
+           {
+             q2=points[p]+ht[i];
+             q1=q2+ht[i];
+             q3=points[p]-ht[i];
+             q4=q3-ht[i];
+             for (unsigned int comp=0; comp<n_components; ++comp)
+               gradients[p][comp][i]=(-  value(q1, comp)
+                                      +8*value(q2, comp)
+                                      -8*value(q3, comp)
+                                      +  value(q4, comp))/(12*h);
+           }
+       break;
+      }
+      default:
+           Assert(false, ExcInvalidFormula());
+    }
+}
+
+
+template <int dim>
+AutoDerivativeFunction<dim>::DifferenceFormula
+AutoDerivativeFunction<dim>::get_formula_of_order(const unsigned int ord)
+{
+  switch (ord)
+    {
+      case 0:
+      case 1: return UpwindEuler;
+      case 2: return Euler;
+      case 3:
+      case 4: return FourthOrder;
+      default:
+           Assert(false, ExcNotImplemented());
+    }
+  return Euler;
+}
+
+
+template class AutoDerivativeFunction<1>;
+template class AutoDerivativeFunction<2>;
+template class AutoDerivativeFunction<3>;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.