--- /dev/null
+//---------------------------- 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
+
+
--- /dev/null
+//---------------------------- 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
+
+
--- /dev/null
+//---------------------------- 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);
+
+
+
--- /dev/null
+//---------------------------- 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);
+
+
+