/**
- * 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.
+ * 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: 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
+ * 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'}\| }.$
* \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'}\|}
+ * $ 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).$
+ * 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.
+ * 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
* gradients themselves.
*
*
- * @sect2{Refinement indicators based on the gradients}
+ * @sect2{Approximation of higher derivatives}
*
- * If you would like to base a refinement criterion upon this
- * approximation of the gradient, you will have to scale the results
+ * 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
* 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 GradientEstimator
+class DerivativeApproximation
{
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
+ * 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 estimate (const DoFHandler<dim> &dof,
+ static void
+ approximate_gradient (const DoFHandler<dim> &dof,
const Vector<double> &solution,
- Vector<float> &error_per_cell);
+ 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
*/
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
typedef pair<unsigned int,unsigned int> IndexInterval;
/**
- * Compute the error estimator on
- * the cells in the range given
- * by the third parameter.
+ * 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 <int dim>
- static void estimate_threaded (const DoFHandler<dim> &dof,
- const Vector<double> &solution,
- const IndexInterval &index_interval,
- Vector<float> &error_per_cell);
+ 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);
};
/**
- * 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.
+ * 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: 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
+ * 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'}\| }.$
* \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'}\|}
+ * $ 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).$
+ * 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.
+ * 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
* gradients themselves.
*
*
- * @sect2{Refinement indicators based on the gradients}
+ * @sect2{Approximation of higher derivatives}
*
- * If you would like to base a refinement criterion upon this
- * approximation of the gradient, you will have to scale the results
+ * 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
* 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 GradientEstimator
+class DerivativeApproximation
{
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
+ * 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 estimate (const DoFHandler<dim> &dof,
+ static void
+ approximate_gradient (const DoFHandler<dim> &dof,
const Vector<double> &solution,
- Vector<float> &error_per_cell);
+ 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
*/
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
typedef pair<unsigned int,unsigned int> IndexInterval;
/**
- * Compute the error estimator on
- * the cells in the range given
- * by the third parameter.
+ * 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 <int dim>
- static void estimate_threaded (const DoFHandler<dim> &dof,
- const Vector<double> &solution,
- const IndexInterval &index_interval,
- Vector<float> &error_per_cell);
+ 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);
};
#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 <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 <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
+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
-GradientEstimator::estimate (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &error_per_cell)
+DerivativeApproximation::
+approximate_second_derivative (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution,
+ Vector<float> &derivative_norm)
{
- Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
- ExcInvalidVectorLength (error_per_cell.size(),
+ 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());
-#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(),
Threads::ThreadManager 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));
+ Threads::encapsulate
+ (&DerivativeApproximation::
+ template approximate<DerivativeDescription,dim>)
+ .collect_args (dof_handler, solution,
+ index_intervals[i],
+ derivative_norm));
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>
+template <class DerivativeDescription, int dim>
void
-GradientEstimator::estimate_threaded (const DoFHandler<dim> &dof_handler,
+DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution,
const IndexInterval &index_interval,
- Vector<float> &error_per_cell)
+ Vector<float> &derivative_norm)
{
QMidpoint<dim> midpoint_rule;
FEValues<dim> fe_midpoint_value (dof_handler.get_fe(),
midpoint_rule,
- UpdateFlags(update_values |
+ UpdateFlags(DerivativeDescription::update_flags |
update_q_points));
// matrix Y=sum_i y_i y_i^T
// respective entries in the output
// vector:
Vector<float>::iterator
- error_on_this_cell = error_per_cell.begin() + index_interval.first;
+ derivative_norm_on_this_cell
+ = derivative_norm.begin() + index_interval.first;
- DoFHandler<dim>::active_cell_iterator cell, endc;
+ typename DoFHandler<dim>::active_cell_iterator cell, endc;
cell = endc = dof_handler.begin_active();
// (static_cast to avoid warnings
// about unsigned always >=0)
// active neighbors of a cell
// reserve the maximal number of
// active neighbors
- vector<DoFHandler<dim>::active_cell_iterator> 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, ++error_on_this_cell)
+ 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|
- Tensor<1,dim> projected_gradient;
+ // 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
- // 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);
+ // 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
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);
+ typename DoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(n);
if (neighbor->active())
active_neighbors.push_back (neighbor);
else
// present cell
if (dim == 1)
{
- DoFHandler<dim>::cell_iterator neighbor_child = neighbor;
+ typename DoFHandler<dim>::cell_iterator
+ neighbor_child = neighbor;
while (neighbor_child->has_children())
neighbor_child = neighbor_child->child (n==0 ? 1 : 0);
// now loop over all active
// neighbors and collect the
// data we need
- typename vector<DoFHandler<dim>::active_cell_iterator>::const_iterator
+ typename vector<typename 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
+ const typename DoFHandler<dim>::active_cell_iterator
neighbor = *neighbor_ptr;
// reinit fe values object...
// ...and get the value of the
// solution...
- vector<double> neighbor_midpoint_value(1);
- fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+ const typename DerivativeDescription::ProjectedDerivative
+ neighbor_midpoint_value
+ = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
+ solution);
+
// ...and the place where it lives
- Point<dim> neighbor_center = fe_midpoint_value.quadrature_point(0);
+ const Point<dim>
+ neighbor_center = fe_midpoint_value.quadrature_point(0);
// vector for the
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
for (unsigned int j=0; j<dim; ++j)
Y[i][j] += y[i] * y[j];
- // the update the sum
+ // then update the sum
// of difference
// quotients
- projected_gradient += (neighbor_midpoint_value[0] -
- this_midpoint_value[0]) /
- distance *
- y;
+ 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
AssertThrow (determinant(Y) != 0,
ExcInsufficientDirections());
+ // first symmetrize g
+ DerivativeDescription::symmetrize (projected_derivative);
+
// compute Y^-1 g
- Point<dim> gradient;
+ typename DerivativeDescription::Derivative derivative;
Tensor<2,dim> Y_inverse = invert(Y);
- // compute Y^-1 g
- contract (gradient, Y_inverse, projected_gradient);
+ contract (derivative, Y_inverse, projected_derivative);
- *error_on_this_cell = sqrt(gradient.square());
+ *derivative_norm_on_this_cell
+ = DerivativeDescription::derivative_norm (derivative);
};
};
// explicit instantiations
template
void
-GradientEstimator::estimate (const DoFHandler<deal_II_dimension> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &error_per_cell);
+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);
#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 <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 <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
+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
-GradientEstimator::estimate (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &error_per_cell)
+DerivativeApproximation::
+approximate_second_derivative (const DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution,
+ Vector<float> &derivative_norm)
{
- Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
- ExcInvalidVectorLength (error_per_cell.size(),
+ 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());
-#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(),
Threads::ThreadManager 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));
+ Threads::encapsulate
+ (&DerivativeApproximation::
+ template approximate<DerivativeDescription,dim>)
+ .collect_args (dof_handler, solution,
+ index_intervals[i],
+ derivative_norm));
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>
+template <class DerivativeDescription, int dim>
void
-GradientEstimator::estimate_threaded (const DoFHandler<dim> &dof_handler,
+DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution,
const IndexInterval &index_interval,
- Vector<float> &error_per_cell)
+ Vector<float> &derivative_norm)
{
QMidpoint<dim> midpoint_rule;
FEValues<dim> fe_midpoint_value (dof_handler.get_fe(),
midpoint_rule,
- UpdateFlags(update_values |
+ UpdateFlags(DerivativeDescription::update_flags |
update_q_points));
// matrix Y=sum_i y_i y_i^T
// respective entries in the output
// vector:
Vector<float>::iterator
- error_on_this_cell = error_per_cell.begin() + index_interval.first;
+ derivative_norm_on_this_cell
+ = derivative_norm.begin() + index_interval.first;
- DoFHandler<dim>::active_cell_iterator cell, endc;
+ typename DoFHandler<dim>::active_cell_iterator cell, endc;
cell = endc = dof_handler.begin_active();
// (static_cast to avoid warnings
// about unsigned always >=0)
// active neighbors of a cell
// reserve the maximal number of
// active neighbors
- vector<DoFHandler<dim>::active_cell_iterator> 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, ++error_on_this_cell)
+ 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|
- Tensor<1,dim> projected_gradient;
+ // 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
- // 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);
+ // 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
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);
+ typename DoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(n);
if (neighbor->active())
active_neighbors.push_back (neighbor);
else
// present cell
if (dim == 1)
{
- DoFHandler<dim>::cell_iterator neighbor_child = neighbor;
+ typename DoFHandler<dim>::cell_iterator
+ neighbor_child = neighbor;
while (neighbor_child->has_children())
neighbor_child = neighbor_child->child (n==0 ? 1 : 0);
// now loop over all active
// neighbors and collect the
// data we need
- typename vector<DoFHandler<dim>::active_cell_iterator>::const_iterator
+ typename vector<typename 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
+ const typename DoFHandler<dim>::active_cell_iterator
neighbor = *neighbor_ptr;
// reinit fe values object...
// ...and get the value of the
// solution...
- vector<double> neighbor_midpoint_value(1);
- fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+ const typename DerivativeDescription::ProjectedDerivative
+ neighbor_midpoint_value
+ = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
+ solution);
+
// ...and the place where it lives
- Point<dim> neighbor_center = fe_midpoint_value.quadrature_point(0);
+ const Point<dim>
+ neighbor_center = fe_midpoint_value.quadrature_point(0);
// vector for the
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
for (unsigned int j=0; j<dim; ++j)
Y[i][j] += y[i] * y[j];
- // the update the sum
+ // then update the sum
// of difference
// quotients
- projected_gradient += (neighbor_midpoint_value[0] -
- this_midpoint_value[0]) /
- distance *
- y;
+ 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
AssertThrow (determinant(Y) != 0,
ExcInsufficientDirections());
+ // first symmetrize g
+ DerivativeDescription::symmetrize (projected_derivative);
+
// compute Y^-1 g
- Point<dim> gradient;
+ typename DerivativeDescription::Derivative derivative;
Tensor<2,dim> Y_inverse = invert(Y);
- // compute Y^-1 g
- contract (gradient, Y_inverse, projected_gradient);
+ contract (derivative, Y_inverse, projected_derivative);
- *error_on_this_cell = sqrt(gradient.square());
+ *derivative_norm_on_this_cell
+ = DerivativeDescription::derivative_norm (derivative);
};
};
// explicit instantiations
template
void
-GradientEstimator::estimate (const DoFHandler<deal_II_dimension> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &error_per_cell);
+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);
projects / dealii.git / commitdiff