]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make DerivativeApproximation a namespace, rather than a class with all static members.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 8 Feb 2014 13:56:26 +0000 (13:56 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 8 Feb 2014 13:56:26 +0000 (13:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@32438 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/numerics/derivative_approximation.h
deal.II/source/numerics/derivative_approximation.cc
deal.II/source/numerics/derivative_approximation.inst.in

index 0126d3df680828cb4bc395b6da6c487ed92e34e6..0b5acf04edeaccefbe0466e168d1827838bb2dfc 100644 (file)
@@ -124,6 +124,11 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>Changed: DerivativeApproximation used to be a class that only had
+  static members. It is now a namespace.
+  <br>
+  (Wolfgang Bangerth, 2014/02/08)
+
   <li>New: ThreadLocalStorage::clear() clears out all objects allocated on the
   current and all other threads.
   <br>
@@ -172,7 +177,7 @@ inconvenience this causes.
   (Matthias Maier, 2014/02/01)
 
   <li>New/fixed: The ParameterHandler::print_parameters_section
-  method not worked for XML output. There is now a flag
+  method did not work for XML output. There is now a flag
   include_top_level_elements which prints all higher
   subsection elements, default is false.
   For XML output setting this flag to true is required
index cc43eeb4790180d4b475a3b0c165eaad2a343d69..660365f7b1da625b88d986668ac8fb8300c4cc50 100644 (file)
@@ -39,7 +39,7 @@ namespace hp
 
 
 /**
- * This class computes a cell-wise approximation of the norm of a
+ * This namespace provides functions that compute a cell-wise approximation of the norm of a
  * derivative of a finite element field by taking difference quotients
  * between neighboring cells. This is a rather simple but efficient
  * form to get an error indicator, since it can be computed with
@@ -85,7 +85,7 @@ namespace hp
  * @code
  * --------------------------------------------------------
  * An error occurred in line <749> of file <source/numerics/derivative_approximation.cc> in function
- *     static void DerivativeApproximation::approximate(const Mapping<dim,spacedim>&, const DH<dim,spacedim>&, const InputVector&, unsigned int, const
+ *     void DerivativeApproximation::approximate(const Mapping<dim,spacedim>&, const DH<dim,spacedim>&, const InputVector&, unsigned int, const
  *  std::pair<unsigned int, unsigned int>&, Vector<float>&) [with DerivativeDescription = DerivativeApproximation::Gradient<3>, int
  * dim = 3, DH = DoFHandler, InputVector = Vector<double>]
  * The violated condition was:
@@ -158,7 +158,7 @@ namespace hp
  * The formulae for the computation of approximations to the gradient
  * and to the tensor of second derivatives shown above are very much
  * alike. The basic difference is that in one case the finite
- * difference quotiont is a scalar, while in the other case it is a
+ * difference quotient is a scalar, while in the other case it is a
  * vector. For higher derivatives, this would be a tensor of even
  * higher rank. We then have to form the outer product of this
  * difference quotient with the distance vector $y_{KK'}$, symmetrize
@@ -185,9 +185,8 @@ namespace hp
  * @ingroup numerics
  * @author Wolfgang Bangerth, 2000
  */
-class DerivativeApproximation
+namespace DerivativeApproximation
 {
-public:
   /**
    * This function is used to obtain an approximation of the gradient. Pass it
    * the DoF handler object that describes the finite element field, a nodal
@@ -204,7 +203,7 @@ public:
    * locally relevant unknowns.
    */
   template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-  static void
+  void
   approximate_gradient (const Mapping<dim,spacedim>    &mapping,
                         const DH<dim,spacedim>         &dof,
                         const InputVector     &solution,
@@ -216,7 +215,7 @@ public:
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-  static void
+  void
   approximate_gradient (const DH<dim,spacedim>         &dof,
                         const InputVector     &solution,
                         Vector<float>         &derivative_norm,
@@ -240,7 +239,7 @@ public:
    * locally relevant unknowns.
    */
   template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-  static void
+  void
   approximate_second_derivative (const Mapping<dim,spacedim>    &mapping,
                                  const DH<dim,spacedim>         &dof,
                                  const InputVector     &solution,
@@ -252,7 +251,7 @@ public:
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-  static void
+  void
   approximate_second_derivative (const DH<dim,spacedim>         &dof,
                                  const InputVector     &solution,
                                  Vector<float>         &derivative_norm,
@@ -272,7 +271,7 @@ public:
    * locally relevant unknowns.
    */
   template <int dim, template <int, int> class DH, class InputVector, int order, int spacedim>
-  static void
+  void
   approximate_derivative_tensor (const Mapping<dim,spacedim>                           &mapping,
                                  const DH<dim,spacedim>                                &dof,
                                  const InputVector                            &solution,
@@ -284,7 +283,7 @@ public:
    * Same as above, with <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, template <int, int> class DH, class InputVector, int order, int spacedim>
-  static void
+  void
   approximate_derivative_tensor (const DH<dim,spacedim>                                &dof,
                                  const InputVector                            &solution,
                                  const typename DH<dim,spacedim>::active_cell_iterator &cell,
@@ -295,8 +294,8 @@ public:
    * Return the norm of the derivative.
    */
   template <int dim, int order>
-  static double
-  derivative_norm(const Tensor<order,dim> &derivative);
+  double
+  derivative_norm (const Tensor<order,dim> &derivative);
 
   /**
    * Exception
@@ -309,299 +308,8 @@ public:
    * Exception
    */
   DeclException0 (ExcInsufficientDirections);
+}
 
-private:
-
-  /**
-   * The following class is used to describe the data needed to compute the
-   * finite difference approximation to the gradient on a cell. See the
-   * general documentation of this class for more information on
-   * implementational details.
-   *
-   * @author Wolfgang Bangerth, 2000
-   */
-  template <int dim>
-  class Gradient
-  {
-  public:
-    /**
-     * Declare which data fields have to be updated for the function @p
-     * get_projected_derivative to work.
-     */
-    static const UpdateFlags update_flags;
-
-    /**
-     * Declare the data type which holds the derivative described by this
-     * class.
-     */
-    typedef Tensor<1,dim> Derivative;
-
-    /**
-     * Likewise declare the data type that holds the derivative projected to a
-     * certain directions.
-     */
-    typedef double        ProjectedDerivative;
-
-    /**
-     * Given an FEValues object initialized to a cell, and a solution vector,
-     * extract the desired derivative at the first quadrature point (which is
-     * the only one, as we only evaluate the finite element field at the
-     * center of each cell).
-     */
-    template <class InputVector, int spacedim>
-    static ProjectedDerivative
-    get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
-                              const InputVector    &solution,
-                              const unsigned int    component);
-
-    /**
-     * Return the norm of the derivative object. Here, for the gradient, we
-     * choose the Euclidian norm of the gradient vector.
-     */
-    static double derivative_norm (const Derivative &d);
-
-    /**
-     * If for the present derivative order, symmetrization of the derivative
-     * tensor is necessary, then do so on the argument.
-     *
-     * For the first derivatives, no such thing is necessary, so this function
-     * is a no-op.
-     */
-    static void symmetrize (Derivative &derivative_tensor);
-  };
-
-
-
-  /**
-   * The following class is used to describe the data needed to compute the
-   * finite difference approximation to the second derivatives on a cell. See
-   * the general documentation of this class for more information on
-   * implementational details.
-   *
-   * @author Wolfgang Bangerth, 2000
-   */
-  template <int dim>
-  class SecondDerivative
-  {
-  public:
-    /**
-     * Declare which data fields have to be updated for the function @p
-     * get_projected_derivative to work.
-     */
-    static const UpdateFlags update_flags;
-
-    /**
-     * Declare the data type which holds the derivative described by this
-     * class.
-     */
-    typedef Tensor<2,dim> Derivative;
-
-    /**
-     * Likewise declare the data type that holds the derivative projected to a
-     * certain directions.
-     */
-    typedef Tensor<1,dim> ProjectedDerivative;
-
-    /**
-     * Given an FEValues object initialized to a cell, and a solution vector,
-     * extract the desired derivative at the first quadrature point (which is
-     * the only one, as we only evaluate the finite element field at the
-     * center of each cell).
-     */
-    template <class InputVector, int spacedim>
-    static ProjectedDerivative
-    get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
-                              const InputVector    &solution,
-                              const unsigned int    component);
-
-    /**
-     * Return the norm of the derivative object. Here, for the (symmetric)
-     * tensor of second derivatives, we choose the absolute value of the
-     * largest eigenvalue, which is the matrix norm associated to the $l_2$
-     * norm of vectors. It is also the largest value of the curvature of the
-     * solution.
-     */
-    static double derivative_norm (const Derivative &d);
-
-    /**
-     * If for the present derivative order, symmetrization of the derivative
-     * tensor is necessary, then do so on the argument.
-     *
-     * For the second derivatives, each entry of the tensor is set to the mean
-     * of its value and the value of the transpose element.
-     *
-     * Note that this function actually modifies its argument.
-     */
-    static void symmetrize (Derivative &derivative_tensor);
-  };
-
-  template <int dim>
-  class ThirdDerivative
-  {
-  public:
-    /**
-     * Declare which data fields have to be updated for the function @p
-     * get_projected_derivative to work.
-     */
-    static const UpdateFlags update_flags;
-
-    /**
-     * Declare the data type which
-     * holds the derivative described
-     * by this class.
-     */
-    typedef Tensor<3,dim> Derivative;
-
-    /**
-     * Likewise declare the data type that holds the derivative projected to a
-     * certain directions.
-     */
-    typedef Tensor<2,dim> ProjectedDerivative;
-
-    /**
-     * Given an FEValues object initialized to a cell, and a solution vector,
-     * extract the desired derivative at the first quadrature point (which is
-     * the only one, as we only evaluate the finite element field at the
-     * center of each cell).
-     */
-    template <class InputVector, int spacedim>
-    static ProjectedDerivative
-    get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
-                              const InputVector    &solution,
-                              const unsigned int    component);
-
-    /**
-     * Return the norm of the derivative object. Here, for the (symmetric)
-     * tensor of second derivatives, we choose the absolute value of the
-     * largest eigenvalue, which is the matrix norm associated to the $l_2$
-     * norm of vectors. It is also the largest value of the curvature of the
-     * solution.
-     */
-    static double derivative_norm (const Derivative &d);
-
-    /**
-     * If for the present derivative order, symmetrization of the derivative
-     * tensor is necessary, then do so on the argument.
-     *
-     * For the second derivatives, each entry of the tensor is set to the mean
-     * of its value and the value of the transpose element.
-     *
-     * Note that this function actually modifies its argument.
-     */
-    static void symmetrize (Derivative &derivative_tensor);
-  };
-
-  template <int order, int dim>
-  class DerivativeSelector
-  {
-  public:
-    /**
-     * typedef to select the DerivativeDescription corresponding to the
-     * <tt>order</tt>th derivative. In this general template we set an unvalid
-     * typedef to void, the real typedefs have to be specialized.
-     */
-    typedef void DerivDescr;
-
-  };
-
-  template <int dim>
-  class DerivativeSelector<1,dim>
-  {
-  public:
-
-    typedef Gradient<dim> DerivDescr;
-  };
-
-  template <int dim>
-  class DerivativeSelector<2,dim>
-  {
-  public:
-
-    typedef SecondDerivative<dim> DerivDescr;
-  };
-
-  template <int dim>
-  class DerivativeSelector<3,dim>
-  {
-  public:
-
-    typedef ThirdDerivative<dim> DerivDescr;
-  };
-
-
-
-
-private:
-
-  /**
-   * Convenience typedef denoting the range of indices on which a certain
-   * thread shall operate.
-   */
-  typedef std::pair<unsigned int,unsigned int> IndexInterval;
-
-  /**
-   * Kind of the main function of this class. It is called by the public entry
-   * points to this class with the correct template first argument and then
-   * simply calls the @p approximate function, after setting up several
-   * threads and doing some administration that is independent of the actual
-   * derivative to be computed.
-   *
-   * The @p component argument denotes which component of the solution vector
-   * we are to work on.
-   */
-  template <class DerivativeDescription, int dim,
-            template <int, int> class DH, class InputVector, int spacedim>
-  static void
-  approximate_derivative (const Mapping<dim,spacedim>    &mapping,
-                          const DH<dim,spacedim>         &dof,
-                          const InputVector     &solution,
-                          const unsigned int     component,
-                          Vector<float>         &derivative_norm);
-
-  /**
-   * Compute the derivative approximation on a given cell.  Fill the @p
-   * derivative_norm vector with the norm of the computed derivative tensors
-   * on the cell.
-   */
-  template <class DerivativeDescription, int dim,
-            template <int, int> class DH, class InputVector, int spacedim>
-  static void
-  approximate (SynchronousIterators<std_cxx1x::tuple<typename DH<dim,spacedim>::active_cell_iterator,
-               Vector<float>::iterator> > const &cell,
-               const Mapping<dim,spacedim>    &mapping,
-               const DH<dim,spacedim>         &dof,
-               const InputVector     &solution,
-               const unsigned int     component);
-
-  /**
-   * Compute the derivative approximation on one cell. This computes the full
-   * derivative tensor.
-   */
-  template <class DerivativeDescription, int dim,
-            template <int, int> class DH, class InputVector, int spacedim>
-  static void
-  approximate_cell (const Mapping<dim,spacedim>                            &mapping,
-                    const DH<dim,spacedim>                                 &dof,
-                    const InputVector                             &solution,
-                    const unsigned int                             component,
-                    const typename DH<dim,spacedim>::active_cell_iterator  &cell,
-                    typename DerivativeDescription::Derivative    &derivative);
-};
-
-
-/* -------------- declaration of explicit specializations ------------- */
-
-template <>
-double
-DerivativeApproximation::SecondDerivative<1>::derivative_norm (const Derivative &d);
-
-template <>
-double
-DerivativeApproximation::SecondDerivative<2>::derivative_norm (const Derivative &d);
-
-template <>
-double
-DerivativeApproximation::SecondDerivative<3>::derivative_norm (const Derivative &d);
 
 
 DEAL_II_NAMESPACE_CLOSE
index 76ffcdbac437a1b869cf29ec6724d6086a8c8746..75a4578557aa114bd6d1a36dd4cf96d24d3c4d16 100644 (file)
@@ -45,832 +45,1083 @@ DEAL_II_NAMESPACE_OPEN
 
 
 
-template <typename T>
-static inline T sqr (const T t)
+namespace
 {
-  return t*t;
+  template <typename T>
+  inline T sqr (const T t)
+  {
+    return t*t;
+  }
 }
 
+// --------------- First the classes and functions that describe individual
+// --------------- derivatives
 
+namespace DerivativeApproximation
+{
+  namespace internal
+  {
+    /**
+     * The following class is used to describe the data needed to compute the
+     * finite difference approximation to the gradient on a cell. See the
+     * general documentation of this class for more information on
+     * implementational details.
+     *
+     * @author Wolfgang Bangerth, 2000
+     */
+    template <int dim>
+    class Gradient
+    {
+      public:
+        /**
+         * Declare which data fields have to be updated for the function @p
+         * get_projected_derivative to work.
+         */
+        static const UpdateFlags update_flags;
+
+        /**
+         * Declare the data type which holds the derivative described by this
+         * class.
+         */
+        typedef Tensor<1,dim> Derivative;
+
+        /**
+         * Likewise declare the data type that holds the derivative projected to a
+         * certain directions.
+         */
+        typedef double        ProjectedDerivative;
+
+        /**
+         * Given an FEValues object initialized to a cell, and a solution vector,
+         * extract the desired derivative at the first quadrature point (which is
+         * the only one, as we only evaluate the finite element field at the
+         * center of each cell).
+         */
+        template <class InputVector, int spacedim>
+        static ProjectedDerivative
+        get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
+                                  const InputVector    &solution,
+                                  const unsigned int    component);
+
+        /**
+         * Return the norm of the derivative object. Here, for the gradient, we
+         * choose the Euclidian norm of the gradient vector.
+         */
+        static double derivative_norm (const Derivative &d);
+
+        /**
+         * If for the present derivative order, symmetrization of the derivative
+         * tensor is necessary, then do so on the argument.
+         *
+         * For the first derivatives, no such thing is necessary, so this function
+         * is a no-op.
+         */
+        static void symmetrize (Derivative &derivative_tensor);
+    };
 
-// static variables
-template <int dim>
-const UpdateFlags DerivativeApproximation::Gradient<dim>::update_flags = update_values;
-
-template <int dim>
-const UpdateFlags DerivativeApproximation::SecondDerivative<dim>::update_flags = update_gradients;
+    // static variables
+    template <int dim>
+    const UpdateFlags Gradient<dim>::update_flags = update_values;
 
-template <int dim>
-const UpdateFlags DerivativeApproximation::ThirdDerivative<dim>::update_flags = update_hessians;
 
+    template <int dim>
+    template <class InputVector, int spacedim>
+    inline
+    typename Gradient<dim>::ProjectedDerivative
+    Gradient<dim>::
+    get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
+        const InputVector    &solution,
+        const unsigned int    component)
+        {
+        if (fe_values.get_fe().n_components() == 1)
+          {
+            std::vector<ProjectedDerivative> values (1);
+            fe_values.get_function_values (solution, values);
+            return values[0];
+          }
+        else
+          {
+            std::vector<Vector<double> > values
+            (1, Vector<double>(fe_values.get_fe().n_components()));
+            fe_values.get_function_values (solution, values);
+            return values[0](component);
+          }
+        }
 
 
-// Dummy structures and dummy function used for WorkStream
-namespace internal
-{
-  namespace Assembler
-  {
-    struct Scratch
-    {
-      Scratch() {}
-    };
 
-    struct CopyData
+    template <int dim>
+    inline
+    double
+    Gradient<dim>::derivative_norm (const Derivative &d)
     {
-      CopyData() {}
-    };
-  }
-}
+      double s = 0;
+      for (unsigned int i=0; i<dim; ++i)
+        s += d[i]*d[i];
+      return std::sqrt(s);
+    }
 
 
 
-template <int dim>
-template <class InputVector, int spacedim>
-inline
-typename DerivativeApproximation::Gradient<dim>::ProjectedDerivative
-DerivativeApproximation::Gradient<dim>::
-get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
-                          const InputVector    &solution,
-                          const unsigned int    component)
-{
-  if (fe_values.get_fe().n_components() == 1)
+    template <int dim>
+    inline
+    void
+    Gradient<dim>::symmetrize (Derivative &)
     {
-      std::vector<ProjectedDerivative> values (1);
-      fe_values.get_function_values (solution, values);
-      return values[0];
+      // nothing to do here
     }
-  else
+
+
+
+    /**
+     * The following class is used to describe the data needed to compute the
+     * finite difference approximation to the second derivatives on a cell. See
+     * the general documentation of this class for more information on
+     * implementational details.
+     *
+     * @author Wolfgang Bangerth, 2000
+     */
+    template <int dim>
+    class SecondDerivative
     {
-      std::vector<Vector<double> > values
-      (1, Vector<double>(fe_values.get_fe().n_components()));
-      fe_values.get_function_values (solution, values);
-      return values[0](component);
-    }
-}
+      public:
+        /**
+         * Declare which data fields have to be updated for the function @p
+         * get_projected_derivative to work.
+         */
+        static const UpdateFlags update_flags;
+
+        /**
+         * Declare the data type which holds the derivative described by this
+         * class.
+         */
+        typedef Tensor<2,dim> Derivative;
+
+        /**
+         * Likewise declare the data type that holds the derivative projected to a
+         * certain directions.
+         */
+        typedef Tensor<1,dim> ProjectedDerivative;
+
+        /**
+         * Given an FEValues object initialized to a cell, and a solution vector,
+         * extract the desired derivative at the first quadrature point (which is
+         * the only one, as we only evaluate the finite element field at the
+         * center of each cell).
+         */
+        template <class InputVector, int spacedim>
+        static ProjectedDerivative
+        get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
+                                  const InputVector    &solution,
+                                  const unsigned int    component);
+
+        /**
+         * Return the norm of the derivative object. Here, for the (symmetric)
+         * tensor of second derivatives, we choose the absolute value of the
+         * largest eigenvalue, which is the matrix norm associated to the $l_2$
+         * norm of vectors. It is also the largest value of the curvature of the
+         * solution.
+         */
+        static double derivative_norm (const Derivative &d);
+
+        /**
+         * If for the present derivative order, symmetrization of the derivative
+         * tensor is necessary, then do so on the argument.
+         *
+         * For the second derivatives, each entry of the tensor is set to the mean
+         * of its value and the value of the transpose element.
+         *
+         * Note that this function actually modifies its argument.
+         */
+        static void symmetrize (Derivative &derivative_tensor);
+    };
 
+    template <int dim>
+    const UpdateFlags SecondDerivative<dim>::update_flags = update_gradients;
 
 
-template <int dim>
-inline
-double
-DerivativeApproximation::Gradient<dim>::derivative_norm (const Derivative &d)
-{
-  double s = 0;
-  for (unsigned int i=0; i<dim; ++i)
-    s += d[i]*d[i];
-  return std::sqrt(s);
-}
+    template <int dim>
+    template <class InputVector, int spacedim>
+    inline
+    typename SecondDerivative<dim>::ProjectedDerivative
+    SecondDerivative<dim>::
+    get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
+        const InputVector    &solution,
+        const unsigned int    component)
+        {
+        if (fe_values.get_fe().n_components() == 1)
+          {
+            std::vector<ProjectedDerivative> values (1);
+            fe_values.get_function_gradients (solution, values);
+            return values[0];
+          }
+        else
+          {
+            std::vector<std::vector<ProjectedDerivative> > values
+            (1, std::vector<ProjectedDerivative>(fe_values.get_fe().n_components()));
+            fe_values.get_function_gradients (solution, values);
+            return values[0][component];
+          };
+        }
 
 
 
-template <int dim>
-inline
-void
-DerivativeApproximation::Gradient<dim>::symmetrize (Derivative &)
-{
-  // nothing to do here
-}
+    template <>
+    inline
+    double
+    SecondDerivative<1>::
+    derivative_norm (const Derivative &d)
+    {
+      return std::fabs (d[0][0]);
+    }
 
 
 
-template <int dim>
-template <class InputVector, int spacedim>
-inline
-typename DerivativeApproximation::SecondDerivative<dim>::ProjectedDerivative
-DerivativeApproximation::SecondDerivative<dim>::
-get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
-                          const InputVector    &solution,
-                          const unsigned int    component)
-{
-  if (fe_values.get_fe().n_components() == 1)
+    template <>
+    inline
+    double
+    SecondDerivative<2>::
+    derivative_norm (const Derivative &d)
     {
-      std::vector<ProjectedDerivative> values (1);
-      fe_values.get_function_gradients (solution, values);
-      return values[0];
+      // note that d should be a
+      // symmetric 2x2 tensor, so the
+      // eigenvalues are:
+      //
+      // 1/2(a+b\pm\sqrt((a-b)^2+4c^2))
+      //
+      // if the d_11=a, d_22=b,
+      // d_12=d_21=c
+      const double radicand = dealii::sqr(d[0][0] - d[1][1]) +
+          4*dealii::sqr(d[0][1]);
+      const double eigenvalues[2]
+                               = { 0.5*(d[0][0] + d[1][1] + std::sqrt(radicand)),
+                                   0.5*(d[0][0] + d[1][1] - std::sqrt(radicand))
+      };
+
+      return std::max (std::fabs (eigenvalues[0]),
+                       std::fabs (eigenvalues[1]));
     }
-  else
+
+
+
+    template <>
+    inline
+    double
+    SecondDerivative<3>::
+    derivative_norm (const Derivative &d)
     {
-      std::vector<std::vector<ProjectedDerivative> > values
-      (1, std::vector<ProjectedDerivative>(fe_values.get_fe().n_components()));
-      fe_values.get_function_gradients (solution, values);
-      return values[0][component];
-    };
-}
+      /*
+      compute the three eigenvalues of the tensor @p{d} and take the
+      largest. one could use the following maple script to generate C
+      code:
+
+      with(linalg);
+      readlib(C);
+      A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]);
+      E:=eigenvals(A);
+      EE:=vector(3,[E[1],E[2],E[3]]);
+      C(EE);
+
+      Unfortunately, with both optimized and non-optimized output, at some
+      places the code `sqrt(-1.0)' is emitted, and I don't know what
+      Maple intends to do with it. This happens both with Maple4 and
+      Maple5.
+
+      Fortunately, Roger Young provided the following Fortran code, which
+      is transcribed below to C. The code uses an algorithm that uses the
+      invariants of a symmetric matrix. (The translated algorithm is
+      augmented by a test for R>0, since R==0 indicates that all three
+      eigenvalues are equal.)
+
+
+          PROGRAM MAIN
+
+    C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX
+    C (ROGER YOUNG, 2001)
+
+          IMPLICIT NONE
+
+          REAL*8 A11, A12, A13, A22, A23, A33
+          REAL*8 I1, J2, J3, AM
+          REAL*8 S11, S12, S13, S22, S23, S33
+          REAL*8 SS12, SS23, SS13
+          REAL*8 R,R3, XX,YY, THETA
+          REAL*8 A1,A2,A3
+          REAL*8 PI
+          PARAMETER (PI=3.141592653587932384D0)
+          REAL*8 A,B,C, TOL
+          PARAMETER (TOL=1.D-14)
+
+    C DEFINE A TEST MATRIX
+
+          A11 = -1.D0
+          A12 = 5.D0
+          A13 = 3.D0
+          A22 = -2.D0
+          A23 = 0.5D0
+          A33 = 4.D0
+
+
+          I1 = A11 + A22 + A33
+          AM = I1/3.D0
+
+          S11 = A11 - AM
+          S22 = A22 - AM
+          S33 = A33 - AM
+          S12 = A12
+          S13 = A13
+          S23 = A23
+
+          SS12 = S12*S12
+          SS23 = S23*S23
+          SS13 = S13*S13
+
+          J2 = S11*S11 + S22*S22 + S33*S33
+          J2 = J2 + 2.D0*(SS12 + SS23 + SS13)
+          J2 = J2/2.D0
+
+          J3 = S11**3 + S22**3 + S33**3
+          J3 = J3 + 3.D0*S11*(SS12 + SS13)
+          J3 = J3 + 3.D0*S22*(SS12 + SS23)
+          J3 = J3 + 3.D0*S33*(SS13 + SS23)
+          J3 = J3 + 6.D0*S12*S23*S13
+          J3 = J3/3.D0
+
+          R = SQRT(4.D0*J2/3.D0)
+          R3 = R*R*R
+          XX = 4.D0*J3/R3
+
+          YY = 1.D0 - DABS(XX)
+          IF(YY.LE.0.D0)THEN
+             IF(YY.GT.(-TOL))THEN
+                WRITE(6,*)'Equal roots: XX= ',XX
+                A = -(XX/DABS(XX))*SQRT(J2/3.D0)
+                B = AM + A
+                C = AM - 2.D0*A
+                WRITE(6,*)B,' (twice) ',C
+                STOP
+             ELSE
+                WRITE(6,*)'Error: XX= ',XX
+                STOP
+             ENDIF
+          ENDIF
+
+          THETA = (ACOS(XX))/3.D0
+
+          A1 = AM + R*COS(THETA)
+          A2 = AM + R*COS(THETA + 2.D0*PI/3.D0)
+          A3 = AM + R*COS(THETA + 4.D0*PI/3.D0)
+
+          WRITE(6,*)A1,A2,A3
+
+          STOP
+          END
+
+       */
+
+      const double am = trace(d) / 3.;
+
+      // s := d - trace(d) I
+      Tensor<2,3> s = d;
+      for (unsigned int i=0; i<3; ++i)
+        s[i][i] -= am;
+
+      const double ss01 = s[0][1] * s[0][1],
+          ss12 = s[1][2] * s[1][2],
+          ss02 = s[0][2] * s[0][2];
+
+      const double J2 = (s[0][0]*s[0][0] + s[1][1]*s[1][1] + s[2][2]*s[2][2]
+                                                                          + 2 * (ss01 + ss02 + ss12))  / 2.;
+      const double J3 = (std::pow(s[0][0],3) + std::pow(s[1][1],3) + std::pow(s[2][2],3)
+      + 3. * s[0][0] * (ss01 + ss02)
+      + 3. * s[1][1] * (ss01 + ss12)
+      + 3. * s[2][2] * (ss02 + ss12)
+      + 6. * s[0][1] * s[0][2] * s[1][2]) / 3.;
+
+      const double R  = std::sqrt (4. * J2 / 3.);
+
+      double EE[3] = { 0, 0, 0 };
+      // the eigenvalues are away from
+      // @p{am} in the order of R. thus,
+      // if R<<AM, then we have the
+      // degenerate case with three
+      // identical eigenvalues. check
+      // this first
+      if (R <= 1e-14*std::fabs(am))
+        EE[0] = EE[1] = EE[2] = am;
+      else
+        {
+          // at least two eigenvalues are
+          // distinct
+          const double R3 = R*R*R;
+          const double XX = 4. * J3 / R3;
+          const double YY = 1. - std::fabs(XX);
+
+          Assert (YY > -1e-14, ExcInternalError());
+
+          if (YY < 0)
+            {
+              // two roots are equal
+              const double a = (XX>0 ? -1. : 1.) * R / 2;
+              EE[0] = EE[1] = am + a;
+              EE[2] = am - 2.*a;
+            }
+          else
+            {
+              const double theta = std::acos(XX) / 3.;
+              EE[0] = am + R*std::cos(theta);
+              EE[1] = am + R*std::cos(theta + 2./3.*numbers::PI);
+              EE[2] = am + R*std::cos(theta + 4./3.*numbers::PI);
+            };
+        };
 
+      return std::max (std::fabs (EE[0]),
+                       std::max (std::fabs (EE[1]),
+                                 std::fabs (EE[2])));
+    }
 
 
-template <>
-inline
-double
-DerivativeApproximation::SecondDerivative<1>::
-derivative_norm (const Derivative &d)
-{
-  return std::fabs (d[0][0]);
-}
 
+    template <int dim>
+    inline
+    double
+    SecondDerivative<dim>::
+    derivative_norm (const Derivative &)
+    {
+      // computing the spectral norm is
+      // not so simple in general. it is
+      // feasible for dim==3 as shown
+      // above, since then there are
+      // still closed form expressions of
+      // the roots of the characteristic
+      // polynomial, and they can easily
+      // be computed using
+      // maple. however, for higher
+      // dimensions, some other method
+      // needs to be employed. maybe some
+      // steps of the power method would
+      // suffice?
+      Assert (false, ExcNotImplemented());
+      return 0;
+    }
 
 
-template <>
-inline
-double
-DerivativeApproximation::SecondDerivative<2>::
-derivative_norm (const Derivative &d)
-{
-  // note that d should be a
-  // symmetric 2x2 tensor, so the
-  // eigenvalues are:
-  //
-  // 1/2(a+b\pm\sqrt((a-b)^2+4c^2))
-  //
-  // if the d_11=a, d_22=b,
-  // d_12=d_21=c
-  const double radicand = dealii::sqr(d[0][0] - d[1][1]) +
-                          4*dealii::sqr(d[0][1]);
-  const double eigenvalues[2]
-    = { 0.5*(d[0][0] + d[1][1] + std::sqrt(radicand)),
-        0.5*(d[0][0] + d[1][1] - std::sqrt(radicand))
-      };
 
-  return std::max (std::fabs (eigenvalues[0]),
-                   std::fabs (eigenvalues[1]));
-}
+    template <int dim>
+    inline
+    void
+    SecondDerivative<dim>::symmetrize (Derivative &d)
+    {
+      // symmetrize non-diagonal entries
+      for (unsigned int i=0; i<dim; ++i)
+        for (unsigned int j=i+1; j<dim; ++j)
+          {
+            const double s = (d[i][j] + d[j][i]) / 2;
+            d[i][j] = d[j][i] = s;
+          };
+    }
 
 
 
-template <>
-inline
-double
-DerivativeApproximation::SecondDerivative<3>::
-derivative_norm (const Derivative &d)
-{
-  /*
-    compute the three eigenvalues of the tensor @p{d} and take the
-    largest. one could use the following maple script to generate C
-    code:
-
-    with(linalg);
-    readlib(C);
-    A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]);
-    E:=eigenvals(A);
-    EE:=vector(3,[E[1],E[2],E[3]]);
-    C(EE);
-
-    Unfortunately, with both optimized and non-optimized output, at some
-    places the code `sqrt(-1.0)' is emitted, and I don't know what
-    Maple intends to do with it. This happens both with Maple4 and
-    Maple5.
-
-    Fortunately, Roger Young provided the following Fortran code, which
-    is transcribed below to C. The code uses an algorithm that uses the
-    invariants of a symmetric matrix. (The translated algorithm is
-    augmented by a test for R>0, since R==0 indicates that all three
-    eigenvalues are equal.)
-
-
-        PROGRAM MAIN
-
-  C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX
-  C (ROGER YOUNG, 2001)
-
-        IMPLICIT NONE
-
-        REAL*8 A11, A12, A13, A22, A23, A33
-        REAL*8 I1, J2, J3, AM
-        REAL*8 S11, S12, S13, S22, S23, S33
-        REAL*8 SS12, SS23, SS13
-        REAL*8 R,R3, XX,YY, THETA
-        REAL*8 A1,A2,A3
-        REAL*8 PI
-        PARAMETER (PI=3.141592653587932384D0)
-        REAL*8 A,B,C, TOL
-        PARAMETER (TOL=1.D-14)
-
-  C DEFINE A TEST MATRIX
-
-        A11 = -1.D0
-        A12 = 5.D0
-        A13 = 3.D0
-        A22 = -2.D0
-        A23 = 0.5D0
-        A33 = 4.D0
-
-
-        I1 = A11 + A22 + A33
-        AM = I1/3.D0
-
-        S11 = A11 - AM
-        S22 = A22 - AM
-        S33 = A33 - AM
-        S12 = A12
-        S13 = A13
-        S23 = A23
-
-        SS12 = S12*S12
-        SS23 = S23*S23
-        SS13 = S13*S13
-
-        J2 = S11*S11 + S22*S22 + S33*S33
-        J2 = J2 + 2.D0*(SS12 + SS23 + SS13)
-        J2 = J2/2.D0
-
-        J3 = S11**3 + S22**3 + S33**3
-        J3 = J3 + 3.D0*S11*(SS12 + SS13)
-        J3 = J3 + 3.D0*S22*(SS12 + SS23)
-        J3 = J3 + 3.D0*S33*(SS13 + SS23)
-        J3 = J3 + 6.D0*S12*S23*S13
-        J3 = J3/3.D0
-
-        R = SQRT(4.D0*J2/3.D0)
-        R3 = R*R*R
-        XX = 4.D0*J3/R3
-
-        YY = 1.D0 - DABS(XX)
-        IF(YY.LE.0.D0)THEN
-           IF(YY.GT.(-TOL))THEN
-              WRITE(6,*)'Equal roots: XX= ',XX
-              A = -(XX/DABS(XX))*SQRT(J2/3.D0)
-              B = AM + A
-              C = AM - 2.D0*A
-              WRITE(6,*)B,' (twice) ',C
-              STOP
-           ELSE
-              WRITE(6,*)'Error: XX= ',XX
-              STOP
-           ENDIF
-        ENDIF
-
-        THETA = (ACOS(XX))/3.D0
-
-        A1 = AM + R*COS(THETA)
-        A2 = AM + R*COS(THETA + 2.D0*PI/3.D0)
-        A3 = AM + R*COS(THETA + 4.D0*PI/3.D0)
-
-        WRITE(6,*)A1,A2,A3
-
-        STOP
-        END
-
-  */
-
-  const double am = trace(d) / 3.;
-
-  // s := d - trace(d) I
-  Tensor<2,3> s = d;
-  for (unsigned int i=0; i<3; ++i)
-    s[i][i] -= am;
-
-  const double ss01 = s[0][1] * s[0][1],
-               ss12 = s[1][2] * s[1][2],
-               ss02 = s[0][2] * s[0][2];
-
-  const double J2 = (s[0][0]*s[0][0] + s[1][1]*s[1][1] + s[2][2]*s[2][2]
-                     + 2 * (ss01 + ss02 + ss12))  / 2.;
-  const double J3 = (std::pow(s[0][0],3) + std::pow(s[1][1],3) + std::pow(s[2][2],3)
-                     + 3. * s[0][0] * (ss01 + ss02)
-                     + 3. * s[1][1] * (ss01 + ss12)
-                     + 3. * s[2][2] * (ss02 + ss12)
-                     + 6. * s[0][1] * s[0][2] * s[1][2]) / 3.;
-
-  const double R  = std::sqrt (4. * J2 / 3.);
-
-  double EE[3] = { 0, 0, 0 };
-  // the eigenvalues are away from
-  // @p{am} in the order of R. thus,
-  // if R<<AM, then we have the
-  // degenerate case with three
-  // identical eigenvalues. check
-  // this first
-  if (R <= 1e-14*std::fabs(am))
-    EE[0] = EE[1] = EE[2] = am;
-  else
+    template <int dim>
+    class ThirdDerivative
     {
-      // at least two eigenvalues are
-      // distinct
-      const double R3 = R*R*R;
-      const double XX = 4. * J3 / R3;
-      const double YY = 1. - std::fabs(XX);
+      public:
+        /**
+         * Declare which data fields have to be updated for the function @p
+         * get_projected_derivative to work.
+         */
+        static const UpdateFlags update_flags;
+
+        /**
+         * Declare the data type which
+         * holds the derivative described
+         * by this class.
+         */
+        typedef Tensor<3,dim> Derivative;
+
+        /**
+         * Likewise declare the data type that holds the derivative projected to a
+         * certain directions.
+         */
+        typedef Tensor<2,dim> ProjectedDerivative;
+
+        /**
+         * Given an FEValues object initialized to a cell, and a solution vector,
+         * extract the desired derivative at the first quadrature point (which is
+         * the only one, as we only evaluate the finite element field at the
+         * center of each cell).
+         */
+        template <class InputVector, int spacedim>
+        static ProjectedDerivative
+        get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
+                                  const InputVector    &solution,
+                                  const unsigned int    component);
+
+        /**
+         * Return the norm of the derivative object. Here, for the (symmetric)
+         * tensor of second derivatives, we choose the absolute value of the
+         * largest eigenvalue, which is the matrix norm associated to the $l_2$
+         * norm of vectors. It is also the largest value of the curvature of the
+         * solution.
+         */
+        static double derivative_norm (const Derivative &d);
+
+        /**
+         * If for the present derivative order, symmetrization of the derivative
+         * tensor is necessary, then do so on the argument.
+         *
+         * For the second derivatives, each entry of the tensor is set to the mean
+         * of its value and the value of the transpose element.
+         *
+         * Note that this function actually modifies its argument.
+         */
+        static void symmetrize (Derivative &derivative_tensor);
+    };
 
-      Assert (YY > -1e-14, ExcInternalError());
+    template <int dim>
+    const UpdateFlags ThirdDerivative<dim>::update_flags = update_hessians;
 
-      if (YY < 0)
+
+    template <int dim>
+    template <class InputVector, int spacedim>
+    inline
+    typename ThirdDerivative<dim>::ProjectedDerivative
+    ThirdDerivative<dim>::
+    get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
+        const InputVector    &solution,
+        const unsigned int    component)
         {
-          // two roots are equal
-          const double a = (XX>0 ? -1. : 1.) * R / 2;
-          EE[0] = EE[1] = am + a;
-          EE[2] = am - 2.*a;
+        if (fe_values.get_fe().n_components() == 1)
+          {
+            std::vector<ProjectedDerivative> values (1);
+            fe_values.get_function_hessians (solution, values);
+            return values[0];
+          }
+        else
+          {
+            std::vector<std::vector<ProjectedDerivative> > values
+            (1, std::vector<ProjectedDerivative>(fe_values.get_fe().n_components()));
+            fe_values.get_function_hessians (solution, values);
+            return values[0][component];
+          };
         }
-      else
-        {
-          const double theta = std::acos(XX) / 3.;
-          EE[0] = am + R*std::cos(theta);
-          EE[1] = am + R*std::cos(theta + 2./3.*numbers::PI);
-          EE[2] = am + R*std::cos(theta + 4./3.*numbers::PI);
-        };
-    };
-
-  return std::max (std::fabs (EE[0]),
-                   std::max (std::fabs (EE[1]),
-                             std::fabs (EE[2])));
-}
 
 
 
-template <int dim>
-inline
-double
-DerivativeApproximation::SecondDerivative<dim>::
-derivative_norm (const Derivative &)
-{
-  // computing the spectral norm is
-  // not so simple in general. it is
-  // feasible for dim==3 as shown
-  // above, since then there are
-  // still closed form expressions of
-  // the roots of the characteristic
-  // polynomial, and they can easily
-  // be computed using
-  // maple. however, for higher
-  // dimensions, some other method
-  // needs to be employed. maybe some
-  // steps of the power method would
-  // suffice?
-  Assert (false, ExcNotImplemented());
-  return 0;
-}
+    template <>
+    inline
+    double
+    ThirdDerivative<1>::
+    derivative_norm (const Derivative &d)
+    {
+      return std::fabs (d[0][0][0]);
+    }
 
 
 
-template <int dim>
-inline
-void
-DerivativeApproximation::SecondDerivative<dim>::symmetrize (Derivative &d)
-{
-  // symmetrize non-diagonal entries
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=i+1; j<dim; ++j)
-      {
-        const double s = (d[i][j] + d[j][i]) / 2;
-        d[i][j] = d[j][i] = s;
-      };
-}
+    template <int dim>
+    inline
+    double
+    ThirdDerivative<dim>::
+    derivative_norm (const Derivative &d)
+    {
+      // return the Frobenius-norm. this is a
+      // member function of Tensor<rank_,dim>
+      return d.norm();
+    }
 
 
-template <int dim>
-template <class InputVector, int spacedim>
-inline
-typename DerivativeApproximation::ThirdDerivative<dim>::ProjectedDerivative
-DerivativeApproximation::ThirdDerivative<dim>::
-get_projected_derivative (const FEValues<dim,spacedim>  &fe_values,
-                          const InputVector    &solution,
-                          const unsigned int    component)
-{
-  if (fe_values.get_fe().n_components() == 1)
+    template <int dim>
+    inline
+    void
+    ThirdDerivative<dim>::symmetrize (Derivative &d)
     {
-      std::vector<ProjectedDerivative> values (1);
-      fe_values.get_function_hessians (solution, values);
-      return values[0];
+      // symmetrize non-diagonal entries
+
+      // first do it in the case, that i,j,k are
+      // pairwise different (which can onlky happen
+      // in dim >= 3)
+      for (unsigned int i=0; i<dim; ++i)
+        for (unsigned int j=i+1; j<dim; ++j)
+          for (unsigned int k=j+1; k<dim; ++k)
+            {
+              const double s = (d[i][j][k] +
+                  d[i][k][j] +
+                  d[j][i][k] +
+                  d[j][k][i] +
+                  d[k][i][j] +
+                  d[k][j][i]) / 6;
+              d[i][j][k]
+                      = d[i][k][j]
+                                = d[j][i][k]
+                                          = d[j][k][i]
+                                                    = d[k][i][j]
+                                                              = d[k][j][i]
+                                                                        = s;
+            }
+      // now do the case, where two indices are
+      // equal
+      for (unsigned int i=0; i<dim; ++i)
+        for (unsigned int j=i+1; j<dim; ++j)
+          {
+            // case 1: index i (lower one) is
+            // double
+            const double s = (d[i][i][j] +
+                d[i][j][i] +
+                d[j][i][i] ) / 3;
+            d[i][i][j]
+                    = d[i][j][i]
+                              = d[j][i][i]
+                                        = s;
+
+            // case 2: index j (higher one) is
+            // double
+            const double t = (d[i][j][j] +
+                d[j][i][j] +
+                d[j][j][i] ) / 3;
+            d[i][j][j]
+                    = d[j][i][j]
+                              = d[j][j][i]
+                                        = t;
+          }
     }
-  else
+
+
+    template <int order, int dim>
+    class DerivativeSelector
     {
-      std::vector<std::vector<ProjectedDerivative> > values
-      (1, std::vector<ProjectedDerivative>(fe_values.get_fe().n_components()));
-      fe_values.get_function_hessians (solution, values);
-      return values[0][component];
+      public:
+        /**
+         * typedef to select the DerivativeDescription corresponding to the
+         * <tt>order</tt>th derivative. In this general template we set an unvalid
+         * typedef to void, the real typedefs have to be specialized.
+         */
+        typedef void DerivDescr;
+
     };
-}
 
+    template <int dim>
+    class DerivativeSelector<1,dim>
+    {
+      public:
 
+        typedef Gradient<dim> DerivDescr;
+    };
 
-template <>
-inline
-double
-DerivativeApproximation::ThirdDerivative<1>::
-derivative_norm (const Derivative &d)
-{
-  return std::fabs (d[0][0][0]);
-}
+    template <int dim>
+    class DerivativeSelector<2,dim>
+    {
+      public:
 
+        typedef SecondDerivative<dim> DerivDescr;
+    };
 
+    template <int dim>
+    class DerivativeSelector<3,dim>
+    {
+      public:
 
-template <int dim>
-inline
-double
-DerivativeApproximation::ThirdDerivative<dim>::
-derivative_norm (const Derivative &d)
-{
-  // return the Frobenius-norm. this is a
-  // member function of Tensor<rank_,dim>
-  return d.norm();
+        typedef ThirdDerivative<dim> DerivDescr;
+    };
+  }
 }
 
-
-template <int dim>
-inline
-void
-DerivativeApproximation::ThirdDerivative<dim>::symmetrize (Derivative &d)
+// Dummy structures and dummy function used for WorkStream
+namespace DerivativeApproximation
 {
-  // symmetrize non-diagonal entries
-
-  // first do it in the case, that i,j,k are
-  // pairwise different (which can onlky happen
-  // in dim >= 3)
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=i+1; j<dim; ++j)
-      for (unsigned int k=j+1; k<dim; ++k)
-        {
-          const double s = (d[i][j][k] +
-                            d[i][k][j] +
-                            d[j][i][k] +
-                            d[j][k][i] +
-                            d[k][i][j] +
-                            d[k][j][i]) / 6;
-          d[i][j][k]
-            = d[i][k][j]
-              = d[j][i][k]
-                = d[j][k][i]
-                  = d[k][i][j]
-                    = d[k][j][i]
-                      = s;
-        };
-  // now do the case, where two indices are
-  // equal
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=i+1; j<dim; ++j)
+  namespace internal
+  {
+    namespace Assembler
+    {
+      struct Scratch
       {
-        // case 1: index i (lower one) is
-        // double
-        const double s = (d[i][i][j] +
-                          d[i][j][i] +
-                          d[j][i][i] ) / 3;
-        d[i][i][j]
-          = d[i][j][i]
-            = d[j][i][i]
-              = s;
-
-        // case 2: index j (higher one) is
-        // double
-        const double t = (d[i][j][j] +
-                          d[j][i][j] +
-                          d[j][j][i] ) / 3;
-        d[i][j][j]
-          = d[j][i][j]
-            = d[j][j][i]
-              = t;
+          Scratch() {}
       };
 
+      struct CopyData
+      {
+          CopyData() {}
+      };
+    }
+  }
 }
 
+// ------------------------------- now for the functions that do the
+// ------------------------------- actual work
 
-
-template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::
-approximate_gradient (const Mapping<dim,spacedim>    &mapping,
-                      const DH<dim,spacedim>         &dof_handler,
-                      const InputVector     &solution,
-                      Vector<float>         &derivative_norm,
-                      const unsigned int     component)
+namespace DerivativeApproximation
 {
-  approximate_derivative<Gradient<dim>,dim> (mapping,
-                                             dof_handler,
-                                             solution,
-                                             component,
-                                             derivative_norm);
-}
+  namespace internal
+  {
+    /**
+     * Compute the derivative approximation on one cell. This computes the full
+     * derivative tensor.
+     */
+    template <class DerivativeDescription, int dim,
+    template <int, int> class DH, class InputVector, int spacedim>
+    void
+    approximate_cell (const Mapping<dim,spacedim>                   &mapping,
+                      const DH<dim,spacedim>                        &dof_handler,
+                      const InputVector                             &solution,
+                      const unsigned int                             component,
+                      const typename DH<dim,spacedim>::active_cell_iterator  &cell,
+                      typename DerivativeDescription::Derivative    &derivative)
+    {
+        QMidpoint<dim> midpoint_rule;
+
+        // create collection objects from
+        // single quadratures, mappings,
+        // and finite elements. if we have
+        // an hp DoFHandler,
+        // dof_handler.get_fe() returns a
+        // collection of which we do a
+        // shallow copy instead
+        const hp::QCollection<dim>       q_collection (midpoint_rule);
+        const hp::FECollection<dim>      fe_collection(dof_handler.get_fe());
+        const hp::MappingCollection<dim> mapping_collection (mapping);
+
+        hp::FEValues<dim> x_fe_midpoint_value (mapping_collection, fe_collection,
+            q_collection,
+            DerivativeDescription::update_flags |
+            update_quadrature_points);
+
+        // matrix Y=sum_i y_i y_i^T
+            Tensor<2,dim> Y;
+
+
+        // vector to hold iterators to all
+            // active neighbors of a cell
+            // reserve the maximal number of
+            // active neighbors
+            std::vector<typename DH<dim,spacedim>::active_cell_iterator> active_neighbors;
+            active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
+                GeometryInfo<dim>::max_children_per_face);
+
+            // vector
+            // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
+            // or related type for higher
+            // derivatives
+            typename DerivativeDescription::Derivative projected_derivative;
+
+            // reinit fe values object...
+            x_fe_midpoint_value.reinit (cell);
+            const FEValues<dim> &fe_midpoint_value
+            = x_fe_midpoint_value.get_present_fe_values();
+
+            // ...and get the value of the
+            // projected derivative...
+            const typename DerivativeDescription::ProjectedDerivative
+            this_midpoint_value
+            = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
+                                                               solution,
+                                                               component);
+            // ...and the place where it lives
+            const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
+
+            // loop over all neighbors and
+            // accumulate the difference
+            // quotients from them. note
+            // that things get a bit more
+            // complicated if the neighbor
+            // is more refined than the
+            // present one
+            //
+            // to make processing simpler,
+            // first collect all neighbor
+            // cells in a vector, and then
+            // collect the data from them
+            GridTools::get_active_neighbors<DH<dim,spacedim> >(cell, active_neighbors);
+
+            // now loop over all active
+            // neighbors and collect the
+            // data we need
+            typename std::vector<typename DH<dim,spacedim>::active_cell_iterator>::const_iterator
+            neighbor_ptr = active_neighbors.begin();
+            for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
+              {
+                const typename DH<dim,spacedim>::active_cell_iterator
+                neighbor = *neighbor_ptr;
+
+                // reinit fe values object...
+                x_fe_midpoint_value.reinit (neighbor);
+                const FEValues<dim> &fe_midpoint_value
+                = x_fe_midpoint_value.get_present_fe_values();
+
+                // ...and get the value of the
+                // solution...
+                const typename DerivativeDescription::ProjectedDerivative
+                neighbor_midpoint_value
+                = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
+                                                                   solution, component);
+
+                // ...and the place where it lives
+                const Point<dim>
+                neighbor_center = fe_midpoint_value.quadrature_point(0);
+
+
+                // vector for the
+                // normalized
+                // direction between
+                // the centers of two
+                // cells
+                Point<dim>   y        = neighbor_center - this_center;
+                const double distance = std::sqrt(y.square());
+                // normalize y
+                y /= distance;
+                // *** note that unlike in
+                // the docs, y denotes the
+                // normalized vector
+                // connecting the centers
+                // of the two cells, rather
+                // than the normal
+                // difference! ***
+
+                // add up the
+                // contribution of
+                // this cell to Y
+                for (unsigned int i=0; i<dim; ++i)
+                  for (unsigned int j=0; j<dim; ++j)
+                    Y[i][j] += y[i] * y[j];
+
+                // then update the sum
+                // of difference
+                // quotients
+                typename DerivativeDescription::ProjectedDerivative
+                projected_finite_difference
+                = (neighbor_midpoint_value -
+                    this_midpoint_value);
+                projected_finite_difference /= distance;
+
+                typename DerivativeDescription::Derivative projected_derivative_update;
+                outer_product (projected_derivative_update,
+                               y,
+                               projected_finite_difference);
+                projected_derivative += projected_derivative_update;
+              };
+
+            // can we determine an
+            // approximation of the
+            // gradient for the present
+            // cell? if so, then we need to
+            // have passed over vectors y_i
+            // which span the whole space,
+            // otherwise we would not have
+            // all components of the
+            // gradient
+            AssertThrow (determinant(Y) != 0,
+                         ExcInsufficientDirections());
+
+            // compute Y^-1 g
+            const Tensor<2,dim> Y_inverse = invert(Y);
+
+            contract (derivative, Y_inverse, projected_derivative);
+
+            // finally symmetrize the derivative
+            DerivativeDescription::symmetrize (derivative);
+    }
 
 
-template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::
-approximate_gradient (const DH<dim,spacedim>         &dof_handler,
-                      const InputVector     &solution,
-                      Vector<float>         &derivative_norm,
-                      const unsigned int     component)
-{
-  approximate_derivative<Gradient<dim>,dim> (StaticMappingQ1<dim>::mapping,
-                                             dof_handler,
-                                             solution,
-                                             component,
-                                             derivative_norm);
-}
 
+    /**
+     * Compute the derivative approximation on a given cell.  Fill the @p
+     * derivative_norm vector with the norm of the computed derivative tensors
+     * on the cell.
+     */
+    template <class DerivativeDescription, int dim,
+    template <int, int> class DH, class InputVector, int spacedim>
+    void
+    approximate (SynchronousIterators<std_cxx1x::tuple<typename DH<dim,spacedim>::active_cell_iterator,Vector<float>::iterator> > const &cell,
+                                          const Mapping<dim,spacedim>                  &mapping,
+                                          const DH<dim,spacedim>                       &dof_handler,
+                                          const InputVector                            &solution,
+                                          const unsigned int                            component)
+    {
+        // if the cell is not locally owned, then there is nothing to do
+        if (std_cxx1x::get<0>(cell.iterators)->is_locally_owned() == false)
+          *std_cxx1x::get<1>(cell.iterators) = 0;
+        else
+          {
+            typename DerivativeDescription::Derivative derivative;
+            // call the function doing the actual
+            // work on this cell
+            approximate_cell<DerivativeDescription,dim,DH,InputVector>
+            (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative);
+            // evaluate the norm and fill the vector
+            //*derivative_norm_on_this_cell
+            *std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative);
+          }
+    }
 
-template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::
-approximate_second_derivative (const Mapping<dim,spacedim>    &mapping,
-                               const DH<dim,spacedim>         &dof_handler,
-                               const InputVector     &solution,
-                               Vector<float>         &derivative_norm,
-                               const unsigned int     component)
-{
-  approximate_derivative<SecondDerivative<dim>,dim> (mapping,
-                                                     dof_handler,
-                                                     solution,
-                                                     component,
-                                                     derivative_norm);
-}
 
+    /**
+     * Kind of the main function of this class. It is called by the public entry
+     * points to this class with the correct template first argument and then
+     * simply calls the @p approximate function, after setting up several
+     * threads and doing some administration that is independent of the actual
+     * derivative to be computed.
+     *
+     * The @p component argument denotes which component of the solution vector
+     * we are to work on.
+     */
+    template <class DerivativeDescription, int dim,
+    template <int, int> class DH, class InputVector, int spacedim>
+    void
+    approximate_derivative (const Mapping<dim,spacedim>    &mapping,
+                            const DH<dim,spacedim>         &dof_handler,
+                            const InputVector     &solution,
+                            const unsigned int     component,
+                            Vector<float>         &derivative_norm)
+    {
+        Assert (derivative_norm.size() == dof_handler.get_tria().n_active_cells(),
+            ExcInvalidVectorLength (derivative_norm.size(),
+                dof_handler.get_tria().n_active_cells()));
+        Assert (component < dof_handler.get_fe().n_components(),
+                ExcIndexRange (component, 0, dof_handler.get_fe().n_components()));
+
+        typedef std_cxx1x::tuple<typename DH<dim,spacedim>::active_cell_iterator,Vector<float>::iterator>
+        Iterators;
+        SynchronousIterators<Iterators> begin(Iterators(dof_handler.begin_active(),
+            derivative_norm.begin())),
+                end(Iterators(dof_handler.end(),
+                    derivative_norm.end()));
+
+        // There is no need for a copier because there is no conflict between threads
+        // to write in derivative_norm. Scratch and CopyData are also useless.
+        WorkStream::run(begin,
+                        end,
+                        static_cast<std_cxx1x::function<void (SynchronousIterators<Iterators> const &,
+                            Assembler::Scratch const &, Assembler::CopyData &)> >
+        (std_cxx1x::bind(&approximate<DerivativeDescription,dim,DH,InputVector,spacedim>,
+                         std_cxx1x::_1,
+                         std_cxx1x::cref(mapping),
+                         std_cxx1x::cref(dof_handler),
+                         std_cxx1x::cref(solution),component)),
+                         std_cxx1x::function<void (internal::Assembler::CopyData const &)> (),
+                         internal::Assembler::Scratch (),internal::Assembler::CopyData ());
+    }
 
-template <int dim, template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::
-approximate_second_derivative (const DH<dim,spacedim>         &dof_handler,
-                               const InputVector     &solution,
-                               Vector<float>         &derivative_norm,
-                               const unsigned int     component)
-{
-  approximate_derivative<SecondDerivative<dim>,dim> (StaticMappingQ1<dim>::mapping,
-                                                     dof_handler,
-                                                     solution,
-                                                     component,
-                                                     derivative_norm);
-}
+  } // namespace internal
 
+} // namespace DerivativeApproximation
 
-template <int dim, template <int, int> class DH, class InputVector, int order, int spacedim>
-void
-DerivativeApproximation::
-approximate_derivative_tensor (const Mapping<dim,spacedim>                           &mapping,
-                               const DH<dim,spacedim>                                &dof,
-                               const InputVector                            &solution,
-                               const typename DH<dim,spacedim>::active_cell_iterator &cell,
-                               Tensor<order,dim>                            &derivative,
-                               const unsigned int                            component)
+
+// ------------------------ finally for the public interface of this namespace
+
+namespace DerivativeApproximation
 {
-  approximate_cell<typename DerivativeSelector<order,dim>::DerivDescr,dim,DH,InputVector>
-  (mapping,
-   dof,
-   solution,
-   component,
-   cell,
-   derivative);
-}
+  template <int dim, template <int, int> class DH, class InputVector, int spacedim>
+  void
+  approximate_gradient (const Mapping<dim,spacedim>    &mapping,
+                        const DH<dim,spacedim>         &dof_handler,
+                        const InputVector     &solution,
+                        Vector<float>         &derivative_norm,
+                        const unsigned int     component)
+  {
+      internal::approximate_derivative<internal::Gradient<dim>,dim> (mapping,
+          dof_handler,
+          solution,
+          component,
+          derivative_norm);
+  }
 
 
+  template <int dim, template <int, int> class DH, class InputVector, int spacedim>
+  void
+  approximate_gradient (const DH<dim,spacedim>         &dof_handler,
+                        const InputVector     &solution,
+                        Vector<float>         &derivative_norm,
+                        const unsigned int     component)
+  {
+      internal::approximate_derivative<internal::Gradient<dim>,dim> (StaticMappingQ1<dim>::mapping,
+          dof_handler,
+          solution,
+          component,
+          derivative_norm);
+  }
 
-template <int dim, template <int, int> class DH, class InputVector, int order, int spacedim>
-void
-DerivativeApproximation::
-approximate_derivative_tensor (const DH<dim,spacedim>                                &dof,
-                               const InputVector                            &solution,
-                               const typename DH<dim,spacedim>::active_cell_iterator &cell,
-                               Tensor<order,dim>                            &derivative,
-                               const unsigned int                            component)
-{
-  // just call the respective function with Q1 mapping
-  approximate_derivative_tensor<dim,DH,InputVector,order,spacedim>
-  (StaticMappingQ1<dim>::mapping,
-   dof,
-   solution,
-   cell,
-   derivative,
-   component);
-}
 
+  template <int dim, template <int, int> class DH, class InputVector, int spacedim>
+  void
+  approximate_second_derivative (const Mapping<dim,spacedim>    &mapping,
+                                 const DH<dim,spacedim>         &dof_handler,
+                                 const InputVector     &solution,
+                                 Vector<float>         &derivative_norm,
+                                 const unsigned int     component)
+  {
+      internal::approximate_derivative<internal::SecondDerivative<dim>,dim> (mapping,
+          dof_handler,
+          solution,
+          component,
+          derivative_norm);
+  }
 
 
-template <class DerivativeDescription, int dim,
-          template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::
-approximate_derivative (const Mapping<dim,spacedim>    &mapping,
-                        const DH<dim,spacedim>         &dof_handler,
-                        const InputVector     &solution,
-                        const unsigned int     component,
-                        Vector<float>         &derivative_norm)
-{
-  Assert (derivative_norm.size() == dof_handler.get_tria().n_active_cells(),
-          ExcInvalidVectorLength (derivative_norm.size(),
-                                  dof_handler.get_tria().n_active_cells()));
-  Assert (component < dof_handler.get_fe().n_components(),
-          ExcIndexRange (component, 0, dof_handler.get_fe().n_components()));
-
-  typedef std_cxx1x::tuple<typename DH<dim,spacedim>::active_cell_iterator,Vector<float>::iterator>
-  Iterators;
-  SynchronousIterators<Iterators> begin(Iterators(dof_handler.begin_active(),
-                                                  derivative_norm.begin())),
-                                  end(Iterators(dof_handler.end(),
-                                                derivative_norm.end()));
-
-  // There is no need for a copier because there is no conflict between threads
-  // to write in derivative_norm. Scratch and CopyData are also useless.
-  WorkStream::run(begin,
-                  end,
-                  static_cast<std_cxx1x::function<void (SynchronousIterators<Iterators> const &,
-                                                        internal::Assembler::Scratch const &,internal::Assembler::CopyData &)> >
-                  (std_cxx1x::bind(&DerivativeApproximation::template approximate<DerivativeDescription,dim,DH,
-                                   InputVector,spacedim>,
-                                   std_cxx1x::_1,
-                                   std_cxx1x::cref(mapping),
-                                   std_cxx1x::cref(dof_handler),
-                                   std_cxx1x::cref(solution),component)),
-                  std_cxx1x::function<void (internal::Assembler::CopyData const &)> (),
-                  internal::Assembler::Scratch (),internal::Assembler::CopyData ());
-}
+  template <int dim, template <int, int> class DH, class InputVector, int spacedim>
+  void
+  approximate_second_derivative (const DH<dim,spacedim>         &dof_handler,
+                                 const InputVector     &solution,
+                                 Vector<float>         &derivative_norm,
+                                 const unsigned int     component)
+  {
+      internal::approximate_derivative<internal::SecondDerivative<dim>,dim> (StaticMappingQ1<dim>::mapping,
+          dof_handler,
+          solution,
+          component,
+          derivative_norm);
+  }
 
 
+  template <int dim, template <int, int> class DH, class InputVector, int order, int spacedim>
+  void
+  approximate_derivative_tensor (const Mapping<dim,spacedim>                           &mapping,
+                                 const DH<dim,spacedim>                                &dof,
+                                 const InputVector                            &solution,
+                                 const typename DH<dim,spacedim>::active_cell_iterator &cell,
+                                 Tensor<order,dim>                            &derivative,
+                                 const unsigned int                            component)
+  {
+      internal::approximate_cell<typename internal::DerivativeSelector<order,dim>::DerivDescr,dim,DH,InputVector>
+      (mapping,
+          dof,
+          solution,
+          component,
+          cell,
+          derivative);
+  }
 
-template <class DerivativeDescription, int dim,
-          template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::approximate (SynchronousIterators<std_cxx1x::tuple<typename DH<dim,spacedim>::active_cell_iterator,Vector<float>::iterator> > const &cell,
-                                      const Mapping<dim,spacedim>                  &mapping,
-                                      const DH<dim,spacedim>                       &dof_handler,
-                                      const InputVector                            &solution,
-                                      const unsigned int                            component)
-{
-  // if the cell is not locally owned, then there is nothing to do
-  if (std_cxx1x::get<0>(cell.iterators)->is_locally_owned() == false)
-    *std_cxx1x::get<1>(cell.iterators) = 0;
-  else
-    {
-      typename DerivativeDescription::Derivative derivative;
-      // call the function doing the actual
-      // work on this cell
-      DerivativeApproximation::template approximate_cell<DerivativeDescription,dim,DH,InputVector>
-      (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative);
-      // evaluate the norm and fill the vector
-      //*derivative_norm_on_this_cell
-      *std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative);
-    }
-}
 
 
+  template <int dim, template <int, int> class DH, class InputVector, int order, int spacedim>
+  void
+  approximate_derivative_tensor (const DH<dim,spacedim>                                &dof,
+                                 const InputVector                            &solution,
+                                 const typename DH<dim,spacedim>::active_cell_iterator &cell,
+                                 Tensor<order,dim>                            &derivative,
+                                 const unsigned int                            component)
+  {
+      // just call the respective function with Q1 mapping
+      approximate_derivative_tensor<dim,DH,InputVector,order,spacedim>
+      (StaticMappingQ1<dim>::mapping,
+          dof,
+          solution,
+          cell,
+          derivative,
+          component);
+  }
 
-template <class DerivativeDescription, int dim,
-          template <int, int> class DH, class InputVector, int spacedim>
-void
-DerivativeApproximation::
-approximate_cell (const Mapping<dim,spacedim>                   &mapping,
-                  const DH<dim,spacedim>                        &dof_handler,
-                  const InputVector                             &solution,
-                  const unsigned int                             component,
-                  const typename DH<dim,spacedim>::active_cell_iterator  &cell,
-                  typename DerivativeDescription::Derivative    &derivative)
-{
-  QMidpoint<dim> midpoint_rule;
-
-  // create collection objects from
-  // single quadratures, mappings,
-  // and finite elements. if we have
-  // an hp DoFHandler,
-  // dof_handler.get_fe() returns a
-  // collection of which we do a
-  // shallow copy instead
-  const hp::QCollection<dim>       q_collection (midpoint_rule);
-  const hp::FECollection<dim>      fe_collection(dof_handler.get_fe());
-  const hp::MappingCollection<dim> mapping_collection (mapping);
-
-  hp::FEValues<dim> x_fe_midpoint_value (mapping_collection, fe_collection,
-                                         q_collection,
-                                         DerivativeDescription::update_flags |
-                                         update_quadrature_points);
-
-  // matrix Y=sum_i y_i y_i^T
-  Tensor<2,dim> Y;
-
-
-  // vector to hold iterators to all
-  // active neighbors of a cell
-  // reserve the maximal number of
-  // active neighbors
-  std::vector<typename DH<dim,spacedim>::active_cell_iterator> active_neighbors;
-  active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
-                            GeometryInfo<dim>::max_children_per_face);
-
-  // vector
-  // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
-  // or related type for higher
-  // derivatives
-  typename DerivativeDescription::Derivative projected_derivative;
-
-  // reinit fe values object...
-  x_fe_midpoint_value.reinit (cell);
-  const FEValues<dim> &fe_midpoint_value
-    = x_fe_midpoint_value.get_present_fe_values();
-
-  // ...and get the value of the
-  // projected derivative...
-  const typename DerivativeDescription::ProjectedDerivative
-  this_midpoint_value
-    = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
-                                                       solution,
-                                                       component);
-  // ...and the place where it lives
-  const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
-
-  // loop over all neighbors and
-  // accumulate the difference
-  // quotients from them. note
-  // that things get a bit more
-  // complicated if the neighbor
-  // is more refined than the
-  // present one
-  //
-  // to make processing simpler,
-  // first collect all neighbor
-  // cells in a vector, and then
-  // collect the data from them
-  GridTools::get_active_neighbors<DH<dim,spacedim> >(cell, active_neighbors);
-
-  // now loop over all active
-  // neighbors and collect the
-  // data we need
-  typename std::vector<typename DH<dim,spacedim>::active_cell_iterator>::const_iterator
-  neighbor_ptr = active_neighbors.begin();
-  for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
-    {
-      const typename DH<dim,spacedim>::active_cell_iterator
-      neighbor = *neighbor_ptr;
-
-      // reinit fe values object...
-      x_fe_midpoint_value.reinit (neighbor);
-      const FEValues<dim> &fe_midpoint_value
-        = x_fe_midpoint_value.get_present_fe_values();
-
-      // ...and get the value of the
-      // solution...
-      const typename DerivativeDescription::ProjectedDerivative
-      neighbor_midpoint_value
-        = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
-                                                           solution, component);
-
-      // ...and the place where it lives
-      const Point<dim>
-      neighbor_center = fe_midpoint_value.quadrature_point(0);
-
-
-      // vector for the
-      // normalized
-      // direction between
-      // the centers of two
-      // cells
-      Point<dim>   y        = neighbor_center - this_center;
-      const double distance = std::sqrt(y.square());
-      // normalize y
-      y /= distance;
-      // *** note that unlike in
-      // the docs, y denotes the
-      // normalized vector
-      // connecting the centers
-      // of the two cells, rather
-      // than the normal
-      // difference! ***
-
-      // add up the
-      // contribution of
-      // this cell to Y
-      for (unsigned int i=0; i<dim; ++i)
-        for (unsigned int j=0; j<dim; ++j)
-          Y[i][j] += y[i] * y[j];
-
-      // then update the sum
-      // of difference
-      // quotients
-      typename DerivativeDescription::ProjectedDerivative
-      projected_finite_difference
-        = (neighbor_midpoint_value -
-           this_midpoint_value);
-      projected_finite_difference /= distance;
-
-      typename DerivativeDescription::Derivative projected_derivative_update;
-      outer_product (projected_derivative_update,
-                     y,
-                     projected_finite_difference);
-      projected_derivative += projected_derivative_update;
-    };
 
-  // can we determine an
-  // approximation of the
-  // gradient for the present
-  // cell? if so, then we need to
-  // have passed over vectors y_i
-  // which span the whole space,
-  // otherwise we would not have
-  // all components of the
-  // gradient
-  AssertThrow (determinant(Y) != 0,
-               ExcInsufficientDirections());
-
-  // compute Y^-1 g
-  const Tensor<2,dim> Y_inverse = invert(Y);
-
-  contract (derivative, Y_inverse, projected_derivative);
-
-  // finally symmetrize the derivative
-  DerivativeDescription::symmetrize (derivative);
-}
 
 
-template <int dim, int order>
-double
-DerivativeApproximation::
-derivative_norm(const Tensor<order,dim> &derivative)
-{
-  return DerivativeSelector<order,dim>::DerivDescr::derivative_norm(derivative);
-}
 
+  template <int dim, int order>
+  double
+  derivative_norm (const Tensor<order,dim> &derivative)
+  {
+    return internal::DerivativeSelector<order,dim>::DerivDescr::derivative_norm(derivative);
+  }
+
+}
 
 
 // --------------------------- explicit instantiations ---------------------
index 6e7ed876abc893734a94beb851560be8a0cf9484..430c8a095f0e7ebdbf738369e7ecc61763218725 100644 (file)
 
 for (deal_II_dimension : DIMENSIONS ; VEC : SERIAL_VECTORS ; DH : DOFHANDLER_TEMPLATES)
 {
+  namespace DerivativeApproximation
+  \{
 template
 void
-DerivativeApproximation::
 approximate_gradient<deal_II_dimension>
 (const Mapping<deal_II_dimension> &mapping,
  const DH<deal_II_dimension> &dof_handler,
@@ -29,7 +30,6 @@ approximate_gradient<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_gradient<deal_II_dimension>
 (const DH<deal_II_dimension> &dof_handler,
  const VEC             &solution,
@@ -38,7 +38,6 @@ approximate_gradient<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_second_derivative<deal_II_dimension>
 (const Mapping<deal_II_dimension> &mapping,
  const DH<deal_II_dimension> &dof_handler,
@@ -48,7 +47,6 @@ approximate_second_derivative<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_second_derivative<deal_II_dimension>
 (const DH<deal_II_dimension> &dof_handler,
  const VEC             &solution,
@@ -57,7 +55,6 @@ approximate_second_derivative<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_derivative_tensor<deal_II_dimension>
 (const Mapping<deal_II_dimension> &    mapping,
  const DH<deal_II_dimension> &dof_handler,
@@ -68,7 +65,6 @@ approximate_derivative_tensor<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_derivative_tensor<deal_II_dimension>
 (const Mapping<deal_II_dimension> &    mapping,
  const DH<deal_II_dimension> &dof_handler,
@@ -79,7 +75,6 @@ approximate_derivative_tensor<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_derivative_tensor<deal_II_dimension>
 (const Mapping<deal_II_dimension> &    mapping,
  const DH<deal_II_dimension> &dof_handler,
@@ -90,7 +85,6 @@ approximate_derivative_tensor<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_derivative_tensor<deal_II_dimension>
 (const DH<deal_II_dimension> &dof_handler,
  const VEC &solution,
@@ -100,7 +94,6 @@ approximate_derivative_tensor<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_derivative_tensor<deal_II_dimension>
 (const DH<deal_II_dimension> &dof_handler,
  const VEC &solution,
@@ -110,7 +103,6 @@ approximate_derivative_tensor<deal_II_dimension>
 
 template
 void
-DerivativeApproximation::
 approximate_derivative_tensor<deal_II_dimension>
 (const DH<deal_II_dimension> &dof_handler,
  const VEC &solution,
@@ -118,40 +110,24 @@ approximate_derivative_tensor<deal_II_dimension>
  Tensor<3,deal_II_dimension> &derivative,
  const unsigned int     component);
 
+  \}
 }
 
 
 for (deal_II_dimension : DIMENSIONS)
 {
+  namespace DerivativeApproximation
+  \{
 template
 double
-DerivativeApproximation::
 derivative_norm(const Tensor<1,deal_II_dimension> &derivative);
 
 template
 double
-DerivativeApproximation::
 derivative_norm(const Tensor<2,deal_II_dimension> &derivative);
 
 template
 double
-DerivativeApproximation::
 derivative_norm(const Tensor<3,deal_II_dimension> &derivative);
-
-
-// static variables
-//
-// on AIX, the linker is unhappy about some missing symbols. they
-// should really be there, but explicitly instantiating them will also
-// not hurt
-template
-const UpdateFlags
-DerivativeApproximation::Gradient<deal_II_dimension>::update_flags;
-
-template
-const UpdateFlags
-DerivativeApproximation::SecondDerivative<deal_II_dimension>::update_flags;
-template
-const UpdateFlags
-DerivativeApproximation::ThirdDerivative<deal_II_dimension>::update_flags;
+  \}
 }

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.