From: Wolfgang Bangerth Date: Fri, 14 Apr 2000 12:24:46 +0000 (+0000) Subject: Move the Gradient estimator from the parameter estimation program to the lib. X-Git-Tag: v8.0.0~20673 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8aab3962221b1e1e75d9959e7e9e1aa757dfdec8;p=dealii.git Move the Gradient estimator from the parameter estimation program to the lib. git-svn-id: https://svn.dealii.org/trunk@2722 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 index 0000000000..7edfc8b742 --- /dev/null +++ b/deal.II/deal.II/include/numerics/derivative_approximation.h @@ -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 +#include +#include + +#include + + + +/** + * 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 + static void estimate (const DoFHandler &dof, + const Vector &solution, + Vector &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 IndexInterval; + + /** + * Compute the error estimator on + * the cells in the range given + * by the third parameter. + */ + template + static void estimate_threaded (const DoFHandler &dof, + const Vector &solution, + const IndexInterval &index_interval, + Vector &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 index 0000000000..7edfc8b742 --- /dev/null +++ b/deal.II/deal.II/include/numerics/gradient_estimator.h @@ -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 +#include +#include + +#include + + + +/** + * 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 + static void estimate (const DoFHandler &dof, + const Vector &solution, + Vector &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 IndexInterval; + + /** + * Compute the error estimator on + * the cells in the range given + * by the third parameter. + */ + template + static void estimate_threaded (const DoFHandler &dof, + const Vector &solution, + const IndexInterval &index_interval, + Vector &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 index 0000000000..c6bd70e131 --- /dev/null +++ b/deal.II/deal.II/source/numerics/derivative_approximation.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_MT +# include +# include +#endif + + + +template +void +GradientEstimator::estimate (const DoFHandler &dof_handler, + const Vector &solution, + Vector &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 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) + .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 +void +GradientEstimator::estimate_threaded (const DoFHandler &dof_handler, + const Vector &solution, + const IndexInterval &index_interval, + Vector &error_per_cell) +{ + QMidpoint midpoint_rule; + FEValues 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::iterator + error_on_this_cell = error_per_cell.begin() + index_interval.first; + + DoFHandler::active_cell_iterator cell, endc; + cell = endc = dof_handler.begin_active(); + // (static_cast to avoid warnings + // about unsigned always >=0) + advance (cell, static_cast(index_interval.first)); + advance (endc, static_cast(index_interval.second)); + + // vector to hold iterators to all + // active neighbors of a cell + // reserve the maximal number of + // active neighbors + vector::active_cell_iterator> active_neighbors; + active_neighbors.reserve (GeometryInfo::faces_per_cell * + GeometryInfo::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 this_midpoint_value(1); + fe_midpoint_value.get_function_values (solution, this_midpoint_value); + // ...and the place where it lives + Point 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::faces_per_cell; ++n) + if (! cell->at_boundary(n)) + { + DoFHandler::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::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::children_per_cell; ++c) + for (unsigned int f=0; f::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::active_cell_iterator>::const_iterator + neighbor_ptr = active_neighbors.begin(); + for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr) + { + const DoFHandler::active_cell_iterator + neighbor = *neighbor_ptr; + + // reinit fe values object... + fe_midpoint_value.reinit (neighbor); + + // ...and get the value of the + // solution... + vector neighbor_midpoint_value(1); + fe_midpoint_value.get_function_values (solution, this_midpoint_value); + // ...and the place where it lives + Point neighbor_center = fe_midpoint_value.quadrature_point(0); + + + // vector for the + // normalized + // direction between + // the centers of two + // cells + Point 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 gradient; + Tensor<2,dim> Y_inverse = invert(Y); + + // compute Y^-1 g + for (unsigned int i=0; i &dof_handler, + const Vector &solution, + Vector &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 index 0000000000..c6bd70e131 --- /dev/null +++ b/deal.II/deal.II/source/numerics/gradient_estimator.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_MT +# include +# include +#endif + + + +template +void +GradientEstimator::estimate (const DoFHandler &dof_handler, + const Vector &solution, + Vector &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 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) + .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 +void +GradientEstimator::estimate_threaded (const DoFHandler &dof_handler, + const Vector &solution, + const IndexInterval &index_interval, + Vector &error_per_cell) +{ + QMidpoint midpoint_rule; + FEValues 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::iterator + error_on_this_cell = error_per_cell.begin() + index_interval.first; + + DoFHandler::active_cell_iterator cell, endc; + cell = endc = dof_handler.begin_active(); + // (static_cast to avoid warnings + // about unsigned always >=0) + advance (cell, static_cast(index_interval.first)); + advance (endc, static_cast(index_interval.second)); + + // vector to hold iterators to all + // active neighbors of a cell + // reserve the maximal number of + // active neighbors + vector::active_cell_iterator> active_neighbors; + active_neighbors.reserve (GeometryInfo::faces_per_cell * + GeometryInfo::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 this_midpoint_value(1); + fe_midpoint_value.get_function_values (solution, this_midpoint_value); + // ...and the place where it lives + Point 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::faces_per_cell; ++n) + if (! cell->at_boundary(n)) + { + DoFHandler::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::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::children_per_cell; ++c) + for (unsigned int f=0; f::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::active_cell_iterator>::const_iterator + neighbor_ptr = active_neighbors.begin(); + for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr) + { + const DoFHandler::active_cell_iterator + neighbor = *neighbor_ptr; + + // reinit fe values object... + fe_midpoint_value.reinit (neighbor); + + // ...and get the value of the + // solution... + vector neighbor_midpoint_value(1); + fe_midpoint_value.get_function_values (solution, this_midpoint_value); + // ...and the place where it lives + Point neighbor_center = fe_midpoint_value.quadrature_point(0); + + + // vector for the + // normalized + // direction between + // the centers of two + // cells + Point 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 gradient; + Tensor<2,dim> Y_inverse = invert(Y); + + // compute Y^-1 g + for (unsigned int i=0; i &dof_handler, + const Vector &solution, + Vector &error_per_cell); + + +