]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Rename these files to derivative_approximation.cc/h.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Jul 2000 08:17:50 +0000 (08:17 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Jul 2000 08:17:50 +0000 (08:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@3198 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/gradient_estimator.h [deleted file]
deal.II/deal.II/source/numerics/gradient_estimator.cc [deleted file]

diff --git a/deal.II/deal.II/include/numerics/gradient_estimator.h b/deal.II/deal.II/include/numerics/gradient_estimator.h
deleted file mode 100644 (file)
index bdff136..0000000
+++ /dev/null
@@ -1,431 +0,0 @@
-//----------------------------  gradient_estimator.h  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2000 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------  gradient_estimator.h  ---------------------------
-#ifndef __deal2__gradient_estimator_h
-#define __deal2__gradient_estimator_h
-
-
-#include <lac/forward_declarations.h>
-#include <grid/forward_declarations.h>
-#include <base/exceptions.h>
-
-#include <utility>
-
-
-
-/**
- * This class computes 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
- * relatively little numerical effort and yet gives a reasonable
- * approximation.
- *
- * The way the difference quotients are computed on cell $K$ is the
- * following (here described for the approximation of the gradient of
- * a finite element field, but see below for higher derivatived): let
- * $K'$ be a neighboring cell, and let $y_{K'}=x_{K'}-x_K$ be the
- * distance vector between the centers of the two cells, then
- *   $ \frac{u_h(x_{K'}) - u_h(x_K)}{ \|y_{K'}\| }$
- * is an approximation of the directional derivative
- *   $ \nabla u(x_K) \cdot \frac{y_{K'}}{ \|y_{K'}\| }.$
- * By multiplying both terms by $\frac{y_{K'}}{ \|y_{K'}\| }$ from the 
- * left and summing over all neighbors $K'$, we obtain
- *   $ \sum_{K'} \left( \frac{y_{K'}}{ \|y_{K'}\|} 
- *                      \frac{y_{K'}^T}{ \|y_{K'}\| } \right) \nabla u(x_K)
- *     \approx
- *     \sum_{K'} \left( \frac{y_{K'}}{ \|y_{K'}\|} 
- *                      \frac{u_h(x_{K'}) - u_h(x_K)}{ \|y_{K'}\| }  \right).$
- *
- * Thus, if the matrix
- *   $ Y =  \sum_{K'} \left( \frac{y_{K'}}{\|y_{K'}\|} 
- *                           \frac{y_{K'}^T}{ \|y_{K'}\| } \right)$ is
- * regular (which is the case when the vectors $y_{K'}$ to all neighbors span
- * the whole space), we can obtain an approximation to the true gradient by
- *   $ \nabla u(x_K)
- *     \approx
- *     Y^{-1} \sum_{K'} \left( \frac{y_{K'}}{\|y_{K'}\|} 
- *                             \frac{u_h(x_{K'}) - u_h(x_K)}{ \|y_{K'}\| }
- *                      \right).$
- * This is a quantity that is easily computed. The value returned for
- * each cell when calling the @p{approximate_gradient} function of
- * this class is the $l_2$ norm of this approximation to the
- * gradient. To make this a useful quantity, you may want to scale
- * each element by the correct power of the respective cell size.
- *
- * The computation of this quantity must fail if a cell has only
- * neighbors for which the direction vectors do not span the whole
- * space. As can easily be verified, this can only happen on very
- * coarse grids, when some cells and all their neighbors have not been
- * refined even once. You should therefore only call the functions of
- * this class if all cells are at least once refined. In practice this
- * is not much of a restriction. If for some cells, the neighbors do
- * not span the whole space, an exception is thrown.
- *
- * Note that for the computation of the quantities of this class, only
- * the values of the finite element field at the centers of the cells
- * are taken. It might therefore only be useful to use this class for
- * discontinuous, piecewise constant elements (i.e. using the
- * @p{FEDG_Q0} class), since all other finite elements can approximate
- * gradients themselves.
- *
- *
- * @sect2{Approximation of higher derivatives}
- *
- * Similar to the reasoning above, approximations to higher
- * derivatives can be computed in a similar fashion. For example, the
- * tensor of second derivatives is approximated by the formula
- *   $ \nabla^2 u(x_K)
- *     \approx
- *     Y^{-1}
- *     \sum_{K'}
- *        \left(
- *           \frac{y_{K'}}{\|y_{K'}\|} \otimes
- *           \frac{\nabla u_h(x_{K'}) - \nabla u_h(x_K)}{ \|y_{K'}\| }
- *        \right),
- *   $ 
- * where $\otimes$ denotes the outer product of two vectors. Note that
- * unlike the true tensor of second derivatives, its approximation is
- * not necessarily symmetric. This is due to the fact that in the
- * derivation, it is not clear whether we shall consider as projected
- * second derivative the term $\nabla^2 u y_{KK'}$ or $y_{KK'}^T
- * \nabla^2 u$. Depending on which choice we take, we obtain one
- * approximation of the tensor of second derivatives or its
- * transpose. To avoid this ambiguity, as result we take the
- * symmetrized form, which is the mean value of the approximation and
- * its transpose.
- *
- * The returned value on each cell is the spectral norm of the
- * approximated tensor of second derivatives, i.e. the largest
- * eigenvalue by absolute value. This equals the largest curvature of
- * the finite element field at each cell, and the spectral norm is the
- * matrix norm associated to the $l_2$ vector norm.
- *
- * Even higher than the second derivative can be obtained along the
- * same lines as exposed above.
- *
- *
- * @sect2{Refinement indicators based on the derivatives}
- *
- * If you would like to base a refinement criterion upon these
- * approximation of the derivatives, you will have to scale the results
- * of this class by an appropriate power of the mesh width. For
- * example, since
- * $\|u-u_h\|^2_{L_2} \le C h^2 \|\nabla u\|^2_{L_2}$, it might be the
- * right thing to scale the indicators as $\eta_K = h \|\nabla u\|_K$,
- * i.e. $\eta_K = h^{1+d/2} \|\nabla u\|_{\infty;K}$, i.e. the right
- * power is $1+d/2$.
- *
- * Likewise, for the second derivative, one should choose a power of
- * the mesh size $h$ one higher than for the gradient.
- *
- *
- * @sect2{Implementation}
- *
- * 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
- * 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
- * it, contract it with the matrix $Y^{-1}$ and compute its norm. To
- * make the implementation simpler and to allow for code reuse, all
- * these operations that are dependent on the actual order of the
- * derivatives to be approximated, as well as the computation of the
- * quantities entering the difference quotient, have been separated
- * into auxiliary nested classes (names @p{Gradient} and
- * @p{SecondDerivative}) and the main algorithm is simply passed one
- * or the other data types and asks them to perform the order
- * dependent operations. The main framework that is independent of
- * this, such as finding all active neighbors, or setting up the
- * matrix $Y$ is done in the main function @p{approximate}.
- *
- * Due to this way of operation, the class may be easily extended for
- * higher oder derivatives than are presently implemented. Basically,
- * only an additional class along the lines of the derivative
- * descriptor classes @p{Gradient} and @p{SecondDerivative} has to be
- * implemented, with the respective typedefs and functions replaced by
- * the appropriate analogues for the derivative that is to be
- * approximated.
- *
- * @author Wolfgang Bangerth, 2000
- */
-class 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 value vector, and
-                                     * receive the cell-wise
-                                     * Euclidian norm of the
-                                     * approximated gradient.
-                                     */
-    template <int dim>
-    static void
-    approximate_gradient (const DoFHandler<dim> &dof,
-                         const Vector<double>  &solution,
-                         Vector<float>         &derivative_norm);
-
-                                    /**
-                                     * This function is the analogue
-                                     * to the one above, computing
-                                     * finite difference
-                                     * approximations of the tensor
-                                     * of second derivatives. Pass it
-                                     * the DoF handler object that
-                                     * describes the finite element
-                                     * field, a nodal value vector,
-                                     * and receive the cell-wise
-                                     * spectral norm of the
-                                     * approximated tensor of second
-                                     * derivatives. The spectral norm
-                                     * is the matrix norm associated
-                                     * to the $l_2$ vector norm.
-                                     */
-    template <int dim>
-    static void
-    approximate_second_derivative (const DoFHandler<dim> &dof,
-                                  const Vector<double>  &solution,
-                                  Vector<float>         &derivative_norm);
-    
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcInvalidVectorLength,
-                   int, int,
-                   << "Vector has length " << arg1 << ", but should have "
-                   << arg2);
-                                    /**
-                                     * 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 = update_values;
-
-                                        /**
-                                         * 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 @p{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).
-                                         */
-       static ProjectedDerivative
-       get_projected_derivative (const FEValues<dim>  &fe_values,
-                                 const Vector<double> &solution);
-    
-                                        /**
-                                         * 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 = update_gradients;
-
-                                        /**
-                                         * 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 @p{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).
-                                         */
-       static ProjectedDerivative
-       get_projected_derivative (const FEValues<dim>  &fe_values,
-                                 const Vector<double> &solution);
-       
-                                        /**
-                                         * 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);
-    };
-    
-                                    /**
-                                     * Convenience typedef denoting
-                                     * the range of indices on which
-                                     * a certain thread shall
-                                     * operate.
-                                     */
-    typedef 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.
-                                     */
-    template <class DerivativeDescription, int dim>
-    static void
-    approximate_derivative (const DoFHandler<dim> &dof,
-                           const Vector<double>  &solution,
-                           Vector<float>         &derivative_norm);
-
-                                    /**
-                                     * Compute the derivative
-                                     * approximation on the cells in
-                                     * the range given by the third
-                                     * parameter.
-                                     */
-    template <class DerivativeDescription, int dim>
-    static void
-    approximate (const DoFHandler<dim> &dof,
-                const Vector<double>  &solution,
-                const IndexInterval   &index_interval,
-                            Vector<float>         &derivative_norm);    
-};
-
-
-#endif
-
-
diff --git a/deal.II/deal.II/source/numerics/gradient_estimator.cc b/deal.II/deal.II/source/numerics/gradient_estimator.cc
deleted file mode 100644 (file)
index 7efe887..0000000
+++ /dev/null
@@ -1,456 +0,0 @@
-//----------------------------  gradient_estimator.cc  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2000 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------  gradient_estimator.cc  ---------------------------
-
-
-#include <base/quadrature_lib.h>
-#include <base/thread_management.h>
-#include <base/multithread_info.h>
-#include <lac/vector.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_handler.h>
-#include <fe/fe.h>
-#include <fe/fe_values.h>
-#include <numerics/gradient_estimator.h>
-
-
-template <typename T>
-static T sqr (const T t)
-{
-  return t*t;
-};
-
-
-
-
-
-template <int dim>
-inline
-typename DerivativeApproximation::Gradient<dim>::ProjectedDerivative
-DerivativeApproximation::Gradient<dim>::
-get_projected_derivative (const FEValues<dim>  &fe_values,
-                         const Vector<double> &solution) 
-{
-  vector<ProjectedDerivative> values (1);
-  fe_values.get_function_values (solution, values);
-  return values[0];
-};
-
-
-
-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 sqrt(s);
-};
-
-
-
-template <int dim>
-inline
-void
-DerivativeApproximation::Gradient<dim>::symmetrize (Derivative &)
-{
-                                  // nothing to do here
-};
-
-
-
-template <int dim>
-inline
-typename DerivativeApproximation::SecondDerivative<dim>::ProjectedDerivative
-DerivativeApproximation::SecondDerivative<dim>::
-get_projected_derivative (const FEValues<dim>  &fe_values,
-                         const Vector<double> &solution) 
-{
-  vector<ProjectedDerivative> values (1);
-  fe_values.get_function_grads (solution, values);
-  return values[0];
-};
-
-
-
-template <>
-inline
-double
-DerivativeApproximation::SecondDerivative<1>::
-derivative_norm (const Derivative &d)
-{
-  return fabs (d[0][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 = sqr(d[0][0] - d[1][1]) + 4*sqr(d[0][1]);
-  const double eigenvalues[2]
-    = { 0.5*(d[0][0] + d[1][1] + sqrt(radicand)),
-       0.5*(d[0][0] + d[1][1] - sqrt(radicand))  };
-  
-  return max (fabs (eigenvalues[0]),
-             fabs (eigenvalues[1]));
-};
-
-
-
-template <int dim>
-inline
-double
-DerivativeApproximation::SecondDerivative<dim>::
-derivative_norm (const Derivative &d)
-{
-                                  // computing the spectral norm is
-                                  // not so simple in general. it is
-                                  // feasible for dim==3, since then
-                                  // there are still closed form
-                                  // expressions of the roots of the
-                                  // third order characteristic
-                                  // polynomial, and they can easily
-                                  // be computed using
-                                  // maple. however, for higher
-                                  // dimensions, some other method
-                                  // needs to be employed.
-  Assert (false, ExcNotImplemented());
-  return 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>
-void 
-DerivativeApproximation::
-approximate_gradient (const DoFHandler<dim> &dof_handler,
-                     const Vector<double>  &solution,
-                     Vector<float>         &derivative_norm)
-{
-  approximate_derivative<Gradient<dim>,dim> (dof_handler,
-                                            solution,
-                                            derivative_norm);
-};
-
-
-
-template <int dim>
-void 
-DerivativeApproximation::
-approximate_second_derivative (const DoFHandler<dim> &dof_handler,
-                              const Vector<double>  &solution,
-                              Vector<float>         &derivative_norm)
-{
-  approximate_derivative<SecondDerivative<dim>,dim> (dof_handler,
-                                                    solution,
-                                                    derivative_norm);
-};
-
-
-
-template <class DerivativeDescription, int dim>
-void 
-DerivativeApproximation::
-approximate_derivative (const DoFHandler<dim> &dof_handler,
-                       const Vector<double>  &solution,
-                       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 (dof_handler.get_fe().n_components() == 1,
-         ExcInternalError());
-
-  const unsigned int n_threads = multithread_info.n_default_threads;
-  vector<IndexInterval> index_intervals
-    = Threads::split_interval (0, dof_handler.get_tria().n_active_cells(),
-                              n_threads);
-  Threads::ThreadManager thread_manager;
-  for (unsigned int i=0; i<n_threads; ++i)
-    Threads::spawn (thread_manager,
-                   Threads::encapsulate
-                   (&DerivativeApproximation::
-                    template approximate<DerivativeDescription,dim>)
-                   .collect_args (dof_handler, solution,
-                                  index_intervals[i],
-                                  derivative_norm));
-  thread_manager.wait ();
-};
-
-
-
-template <class DerivativeDescription, int dim>
-void 
-DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
-                                     const Vector<double>  &solution,
-                                     const IndexInterval   &index_interval,
-                                     Vector<float>         &derivative_norm)
-{
-  QMidpoint<dim> midpoint_rule;
-  FEValues<dim>  fe_midpoint_value (dof_handler.get_fe(),
-                                   midpoint_rule,
-                                   UpdateFlags(DerivativeDescription::update_flags |
-                                               update_q_points));
-  
-                                  // matrix Y=sum_i y_i y_i^T
-  Tensor<2,dim> Y;
-  
-                                  // iterators over all cells and the
-                                  // respective entries in the output
-                                  // vector:
-  Vector<float>::iterator
-    derivative_norm_on_this_cell
-    = derivative_norm.begin() + index_interval.first;
-  
-  typename DoFHandler<dim>::active_cell_iterator cell, endc;
-  cell = endc = dof_handler.begin_active();
-                                  // (static_cast to avoid warnings
-                                  // about unsigned always >=0)
-  advance (cell, static_cast<int>(index_interval.first));
-  advance (endc, static_cast<int>(index_interval.second));
-
-                                  // vector to hold iterators to all
-                                  // active neighbors of a cell
-                                  // reserve the maximal number of
-                                  // active neighbors
-  vector<typename DoFHandler<dim>::active_cell_iterator> active_neighbors;
-  active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
-                           GeometryInfo<dim>::subfaces_per_face);
-
-  for (; cell!=endc; ++cell, ++derivative_norm_on_this_cell)
-    {
-      Y.clear ();
-                                      // 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...
-      fe_midpoint_value.reinit (cell);
-
-                                      // ...and get the value of the
-                                      // projected derivative...
-      const typename DerivativeDescription::ProjectedDerivative
-       this_midpoint_value
-       = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
-                                                   solution);
-                                      // ...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
-      active_neighbors.clear ();
-      for (unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
-       if (! cell->at_boundary(n))
-         {
-           typename DoFHandler<dim>::cell_iterator
-             neighbor = cell->neighbor(n);
-           if (neighbor->active())
-             active_neighbors.push_back (neighbor);
-           else
-             {
-                                                // check children
-                                                // of
-                                                // neighbor. note
-                                                // that in 1d
-                                                // children of
-                                                // the neighbor
-                                                // may be further
-                                                // refined, while
-                                                // they can't in
-                                                // more than one
-                                                // dimension. however,
-                                                // in 1d the case
-                                                // is simpler
-                                                // since we know
-                                                // what children
-                                                // bound to the
-                                                // present cell
-               if (dim == 1)
-                 {
-                   typename DoFHandler<dim>::cell_iterator
-                     neighbor_child = neighbor;
-                   while (neighbor_child->has_children())
-                     neighbor_child = neighbor_child->child (n==0 ? 1 : 0);
-                   
-                   Assert (neighbor_child->neighbor(n==0 ? 1 : 0)==cell,
-                           ExcInternalError());
-                   
-                   active_neighbors.push_back (neighbor_child);
-                 }
-               else
-                                                  // this neighbor has
-                                                  // children. find out
-                                                  // which border to the
-                                                  // present cell
-                 for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-                   for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-                     if (neighbor->child(c)->neighbor(f) == cell)
-                       active_neighbors.push_back (neighbor->child(c));
-             };
-         };
-
-                                      // now loop over all active
-                                      // neighbors and collect the
-                                      // data we need
-      typename vector<typename DoFHandler<dim>::active_cell_iterator>::const_iterator
-       neighbor_ptr = active_neighbors.begin();
-      for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
-       {
-         const typename DoFHandler<dim>::active_cell_iterator
-           neighbor = *neighbor_ptr;
-           
-                                          // reinit fe values object...
-         fe_midpoint_value.reinit (neighbor);
-         
-                                          // ...and get the value of the
-                                          // solution...
-         const typename DerivativeDescription::ProjectedDerivative
-           neighbor_midpoint_value
-           = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
-                                                       solution);
-         
-                                          // ...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 = 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());
-
-                                      // first symmetrize g
-      DerivativeDescription::symmetrize (projected_derivative);
-      
-                                       // compute Y^-1 g
-      typename DerivativeDescription::Derivative derivative;
-      Tensor<2,dim> Y_inverse = invert(Y);
-      
-      contract (derivative, Y_inverse, projected_derivative);
-
-      *derivative_norm_on_this_cell
-       = DerivativeDescription::derivative_norm (derivative);
-    };
-};
-
-
-
-
-// explicit instantiations
-template
-void 
-DerivativeApproximation::
-approximate_gradient (const DoFHandler<deal_II_dimension> &dof_handler,
-                     const Vector<double>  &solution,
-                     Vector<float>         &derivative_norm);
-
-template
-void 
-DerivativeApproximation::
-approximate_second_derivative (const DoFHandler<deal_II_dimension> &dof_handler,
-                              const Vector<double>  &solution,
-                              Vector<float>         &derivative_norm);
-
-
-

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.