]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move the Gradient estimator from the parameter estimation program to the lib.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 14 Apr 2000 12:24:46 +0000 (12:24 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 14 Apr 2000 12:24:46 +0000 (12:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@2722 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/derivative_approximation.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/gradient_estimator.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/derivative_approximation.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/gradient_estimator.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/numerics/derivative_approximation.h b/deal.II/deal.II/include/numerics/derivative_approximation.h
new file mode 100644 (file)
index 0000000..7edfc8b
--- /dev/null
@@ -0,0 +1,147 @@
+//----------------------------  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 estimate of the gradient 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: 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 main 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
+ * #FEDG_Q0# class), since all other finite elements can approximate
+ * gradients themselves.
+ *
+ *
+ * \section{Refinement indicators based on the gradients}
+ *
+ * If you would like to base a refinement criterion upon this
+ * approximation of the gradient, 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 \|\grad u\|^2_{L_2}$, it might be the
+ * right thing to scale the indicators as $\eta_K = h \|\grad u\|_K$,
+ * i.e. $\eta_K = h^{1+d/2} \|\grad u\|_{\infty;K}$, i.e. the right
+ * power is $1+d/2$.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+class GradientEstimator 
+{
+  public:
+                                    /**
+                                     * This is the main function that
+                                     * does what is announced in the
+                                     * general documentation of this
+                                     * class. Pass it the DoF handler
+                                     * object that describes the
+                                     * finite element field, a nodal
+                                     * value vector, and receive the
+                                     * cell-wise norm of the
+                                     * approximated gradient.
+                                     */
+    template <int dim>
+    static void estimate (const DoFHandler<dim> &dof,
+                         const Vector<double>  &solution,
+                         Vector<float>         &error_per_cell);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidVectorLength,
+                   int, int,
+                   << "Vector has length " << arg1 << ", but should have "
+                   << arg2);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInsufficientDirections);
+
+  private:
+                                    /**
+                                     * Convenience typedef denoting
+                                     * the range of indices on which
+                                     * a certain thread shall
+                                     * operate.
+                                     */
+    typedef pair<unsigned int,unsigned int> IndexInterval;
+
+                                    /**
+                                     * Compute the error estimator on
+                                     * the cells in the range given
+                                     * by the third parameter.
+                                     */
+    template <int dim>
+    static void estimate_threaded (const DoFHandler<dim> &dof,
+                                  const Vector<double>  &solution,
+                                  const IndexInterval   &index_interval,
+                                  Vector<float>         &error_per_cell);    
+};
+
+
+#endif
+
+
diff --git a/deal.II/deal.II/include/numerics/gradient_estimator.h b/deal.II/deal.II/include/numerics/gradient_estimator.h
new file mode 100644 (file)
index 0000000..7edfc8b
--- /dev/null
@@ -0,0 +1,147 @@
+//----------------------------  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 estimate of the gradient 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: 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 main 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
+ * #FEDG_Q0# class), since all other finite elements can approximate
+ * gradients themselves.
+ *
+ *
+ * \section{Refinement indicators based on the gradients}
+ *
+ * If you would like to base a refinement criterion upon this
+ * approximation of the gradient, 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 \|\grad u\|^2_{L_2}$, it might be the
+ * right thing to scale the indicators as $\eta_K = h \|\grad u\|_K$,
+ * i.e. $\eta_K = h^{1+d/2} \|\grad u\|_{\infty;K}$, i.e. the right
+ * power is $1+d/2$.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+class GradientEstimator 
+{
+  public:
+                                    /**
+                                     * This is the main function that
+                                     * does what is announced in the
+                                     * general documentation of this
+                                     * class. Pass it the DoF handler
+                                     * object that describes the
+                                     * finite element field, a nodal
+                                     * value vector, and receive the
+                                     * cell-wise norm of the
+                                     * approximated gradient.
+                                     */
+    template <int dim>
+    static void estimate (const DoFHandler<dim> &dof,
+                         const Vector<double>  &solution,
+                         Vector<float>         &error_per_cell);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidVectorLength,
+                   int, int,
+                   << "Vector has length " << arg1 << ", but should have "
+                   << arg2);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInsufficientDirections);
+
+  private:
+                                    /**
+                                     * Convenience typedef denoting
+                                     * the range of indices on which
+                                     * a certain thread shall
+                                     * operate.
+                                     */
+    typedef pair<unsigned int,unsigned int> IndexInterval;
+
+                                    /**
+                                     * Compute the error estimator on
+                                     * the cells in the range given
+                                     * by the third parameter.
+                                     */
+    template <int dim>
+    static void estimate_threaded (const DoFHandler<dim> &dof,
+                                  const Vector<double>  &solution,
+                                  const IndexInterval   &index_interval,
+                                  Vector<float>         &error_per_cell);    
+};
+
+
+#endif
+
+
diff --git a/deal.II/deal.II/source/numerics/derivative_approximation.cc b/deal.II/deal.II/source/numerics/derivative_approximation.cc
new file mode 100644 (file)
index 0000000..c6bd70e
--- /dev/null
@@ -0,0 +1,268 @@
+//----------------------------  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 <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>
+
+#ifdef DEAL_II_USE_MT
+#  include <base/thread_management.h>
+#  include <base/multithread_info.h>
+#endif
+
+
+
+template <int dim>
+void 
+GradientEstimator::estimate (const DoFHandler<dim> &dof_handler,
+                            const Vector<double>  &solution,
+                            Vector<float>         &error_per_cell)
+{
+  Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
+         ExcInvalidVectorLength (error_per_cell.size(),
+                                 dof_handler.get_tria().n_active_cells()));
+  Assert (dof_handler.get_fe().n_components() == 1,
+         ExcInternalError());
+
+#ifdef DEAL_II_USE_MT
+  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);
+  ACE_Thread_Manager thread_manager;
+  for (unsigned int i=0; i<n_threads; ++i)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate (&GradientEstimator::
+                                         template estimate_threaded<dim>)
+                   .collect_args (dof_handler, solution, index_intervals[i],
+                                  error_per_cell));
+  thread_manager.wait ();
+  
+#else
+  estimate_threaded (dof_handler, solution,
+                    make_pair(0U, dof_handler.get_tria().n_active_cells()),
+                    error_per_cell);
+#endif
+};
+
+
+
+template <int dim>
+void 
+GradientEstimator::estimate_threaded (const DoFHandler<dim> &dof_handler,
+                                     const Vector<double>  &solution,
+                                     const IndexInterval   &index_interval,
+                                     Vector<float>         &error_per_cell)
+{
+  QMidpoint<dim> midpoint_rule;
+  FEValues<dim>  fe_midpoint_value (dof_handler.get_fe(),
+                                   midpoint_rule,
+                                   UpdateFlags(update_values |
+                                               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
+    error_on_this_cell = error_per_cell.begin() + index_interval.first;
+  
+  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<DoFHandler<dim>::active_cell_iterator> active_neighbors;
+  active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
+                           GeometryInfo<dim>::subfaces_per_face);
+
+  for (; cell!=endc; ++cell, ++error_on_this_cell)
+    {
+      Y.clear ();
+                                      // vector g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
+      Tensor<1,dim> projected_gradient;
+
+                                      // reinit fe values object...
+      fe_midpoint_value.reinit (cell);
+
+                                      // ...and get the value of the
+                                      // solution...
+      vector<double> this_midpoint_value(1);
+      fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+                                      // ...and the place where it lives
+      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))
+         {
+           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)
+                 {
+                   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
+      vector<DoFHandler<dim>::active_cell_iterator>::const_iterator
+       neighbor_ptr = active_neighbors.begin();
+      for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
+       {
+         const DoFHandler<dim>::active_cell_iterator
+           neighbor = *neighbor_ptr;
+           
+                                          // reinit fe values object...
+         fe_midpoint_value.reinit (neighbor);
+         
+                                          // ...and get the value of the
+                                          // solution...
+         vector<double> neighbor_midpoint_value(1);
+         fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+                                          // ...and the place where it lives
+         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;
+         double     distance = sqrt(y.square());
+                                          // normalize y
+         y /= distance;
+         
+                                          // 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];
+         
+                                          // the update the sum
+                                          // of difference
+                                          // quotients
+         projected_gradient += (neighbor_midpoint_value[0] -
+                                this_midpoint_value[0]) /
+                               distance *
+                               y;
+       };
+
+                                      // 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
+      if (determinant(Y) != 0)
+       {
+                                          // compute Y^-1 g
+         Point<dim> gradient;
+         Tensor<2,dim> Y_inverse = invert(Y);
+         
+                                          // compute Y^-1 g
+         for (unsigned int i=0; i<dim; ++i)
+           for (unsigned int j=0; j<dim; ++j)
+             gradient[i] += Y_inverse[i][j] * projected_gradient[j];
+
+         *error_on_this_cell = sqrt(gradient.square());
+       }
+      else
+                                        // not all search directions
+                                        // available
+       AssertThrow (false, ExcInsufficientDirections());
+    };
+};
+
+
+
+
+// explicit instantiations
+template
+void 
+GradientEstimator::estimate (const DoFHandler<deal_II_dimension> &dof_handler,
+                            const Vector<double>  &solution,
+                            Vector<float>         &error_per_cell);
+
+
+
diff --git a/deal.II/deal.II/source/numerics/gradient_estimator.cc b/deal.II/deal.II/source/numerics/gradient_estimator.cc
new file mode 100644 (file)
index 0000000..c6bd70e
--- /dev/null
@@ -0,0 +1,268 @@
+//----------------------------  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 <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>
+
+#ifdef DEAL_II_USE_MT
+#  include <base/thread_management.h>
+#  include <base/multithread_info.h>
+#endif
+
+
+
+template <int dim>
+void 
+GradientEstimator::estimate (const DoFHandler<dim> &dof_handler,
+                            const Vector<double>  &solution,
+                            Vector<float>         &error_per_cell)
+{
+  Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
+         ExcInvalidVectorLength (error_per_cell.size(),
+                                 dof_handler.get_tria().n_active_cells()));
+  Assert (dof_handler.get_fe().n_components() == 1,
+         ExcInternalError());
+
+#ifdef DEAL_II_USE_MT
+  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);
+  ACE_Thread_Manager thread_manager;
+  for (unsigned int i=0; i<n_threads; ++i)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate (&GradientEstimator::
+                                         template estimate_threaded<dim>)
+                   .collect_args (dof_handler, solution, index_intervals[i],
+                                  error_per_cell));
+  thread_manager.wait ();
+  
+#else
+  estimate_threaded (dof_handler, solution,
+                    make_pair(0U, dof_handler.get_tria().n_active_cells()),
+                    error_per_cell);
+#endif
+};
+
+
+
+template <int dim>
+void 
+GradientEstimator::estimate_threaded (const DoFHandler<dim> &dof_handler,
+                                     const Vector<double>  &solution,
+                                     const IndexInterval   &index_interval,
+                                     Vector<float>         &error_per_cell)
+{
+  QMidpoint<dim> midpoint_rule;
+  FEValues<dim>  fe_midpoint_value (dof_handler.get_fe(),
+                                   midpoint_rule,
+                                   UpdateFlags(update_values |
+                                               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
+    error_on_this_cell = error_per_cell.begin() + index_interval.first;
+  
+  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<DoFHandler<dim>::active_cell_iterator> active_neighbors;
+  active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
+                           GeometryInfo<dim>::subfaces_per_face);
+
+  for (; cell!=endc; ++cell, ++error_on_this_cell)
+    {
+      Y.clear ();
+                                      // vector g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
+      Tensor<1,dim> projected_gradient;
+
+                                      // reinit fe values object...
+      fe_midpoint_value.reinit (cell);
+
+                                      // ...and get the value of the
+                                      // solution...
+      vector<double> this_midpoint_value(1);
+      fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+                                      // ...and the place where it lives
+      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))
+         {
+           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)
+                 {
+                   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
+      vector<DoFHandler<dim>::active_cell_iterator>::const_iterator
+       neighbor_ptr = active_neighbors.begin();
+      for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
+       {
+         const DoFHandler<dim>::active_cell_iterator
+           neighbor = *neighbor_ptr;
+           
+                                          // reinit fe values object...
+         fe_midpoint_value.reinit (neighbor);
+         
+                                          // ...and get the value of the
+                                          // solution...
+         vector<double> neighbor_midpoint_value(1);
+         fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+                                          // ...and the place where it lives
+         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;
+         double     distance = sqrt(y.square());
+                                          // normalize y
+         y /= distance;
+         
+                                          // 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];
+         
+                                          // the update the sum
+                                          // of difference
+                                          // quotients
+         projected_gradient += (neighbor_midpoint_value[0] -
+                                this_midpoint_value[0]) /
+                               distance *
+                               y;
+       };
+
+                                      // 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
+      if (determinant(Y) != 0)
+       {
+                                          // compute Y^-1 g
+         Point<dim> gradient;
+         Tensor<2,dim> Y_inverse = invert(Y);
+         
+                                          // compute Y^-1 g
+         for (unsigned int i=0; i<dim; ++i)
+           for (unsigned int j=0; j<dim; ++j)
+             gradient[i] += Y_inverse[i][j] * projected_gradient[j];
+
+         *error_on_this_cell = sqrt(gradient.square());
+       }
+      else
+                                        // not all search directions
+                                        // available
+       AssertThrow (false, ExcInsufficientDirections());
+    };
+};
+
+
+
+
+// explicit instantiations
+template
+void 
+GradientEstimator::estimate (const DoFHandler<deal_II_dimension> &dof_handler,
+                            const Vector<double>  &solution,
+                            Vector<float>         &error_per_cell);
+
+
+

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.