+++ /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 approximation of the norm of a
- * derivative of a finite element field by taking difference quotients
- * between neighboring cells. This is a rather simple but efficient
- * form to get an error indicator, since it can be computed with
- * relatively little numerical effort and yet gives a reasonable
- * approximation.
- *
- * The way the difference quotients are computed on cell $K$ is the
- * following (here described for the approximation of the gradient of
- * a finite element field, but see below for higher derivatived): let
- * $K'$ be a neighboring cell, and let $y_{K'}=x_{K'}-x_K$ be the
- * distance vector between the centers of the two cells, then
- * $ \frac{u_h(x_{K'}) - u_h(x_K)}{ \|y_{K'}\| }$
- * is an approximation of the directional derivative
- * $ \nabla u(x_K) \cdot \frac{y_{K'}}{ \|y_{K'}\| }.$
- * By multiplying both terms by $\frac{y_{K'}}{ \|y_{K'}\| }$ from the
- * left and summing over all neighbors $K'$, we obtain
- * $ \sum_{K'} \left( \frac{y_{K'}}{ \|y_{K'}\|}
- * \frac{y_{K'}^T}{ \|y_{K'}\| } \right) \nabla u(x_K)
- * \approx
- * \sum_{K'} \left( \frac{y_{K'}}{ \|y_{K'}\|}
- * \frac{u_h(x_{K'}) - u_h(x_K)}{ \|y_{K'}\| } \right).$
- *
- * Thus, if the matrix
- * $ Y = \sum_{K'} \left( \frac{y_{K'}}{\|y_{K'}\|}
- * \frac{y_{K'}^T}{ \|y_{K'}\| } \right)$ is
- * regular (which is the case when the vectors $y_{K'}$ to all neighbors span
- * the whole space), we can obtain an approximation to the true gradient by
- * $ \nabla u(x_K)
- * \approx
- * Y^{-1} \sum_{K'} \left( \frac{y_{K'}}{\|y_{K'}\|}
- * \frac{u_h(x_{K'}) - u_h(x_K)}{ \|y_{K'}\| }
- * \right).$
- * This is a quantity that is easily computed. The value returned for
- * each cell when calling the @p{approximate_gradient} function of
- * this class is the $l_2$ norm of this approximation to the
- * gradient. To make this a useful quantity, you may want to scale
- * each element by the correct power of the respective cell size.
- *
- * The computation of this quantity must fail if a cell has only
- * neighbors for which the direction vectors do not span the whole
- * space. As can easily be verified, this can only happen on very
- * coarse grids, when some cells and all their neighbors have not been
- * refined even once. You should therefore only call the functions of
- * this class if all cells are at least once refined. In practice this
- * is not much of a restriction. If for some cells, the neighbors do
- * not span the whole space, an exception is thrown.
- *
- * Note that for the computation of the quantities of this class, only
- * the values of the finite element field at the centers of the cells
- * are taken. It might therefore only be useful to use this class for
- * discontinuous, piecewise constant elements (i.e. using the
- * @p{FEDG_Q0} class), since all other finite elements can approximate
- * gradients themselves.
- *
- *
- * @sect2{Approximation of higher derivatives}
- *
- * Similar to the reasoning above, approximations to higher
- * derivatives can be computed in a similar fashion. For example, the
- * tensor of second derivatives is approximated by the formula
- * $ \nabla^2 u(x_K)
- * \approx
- * Y^{-1}
- * \sum_{K'}
- * \left(
- * \frac{y_{K'}}{\|y_{K'}\|} \otimes
- * \frac{\nabla u_h(x_{K'}) - \nabla u_h(x_K)}{ \|y_{K'}\| }
- * \right),
- * $
- * where $\otimes$ denotes the outer product of two vectors. Note that
- * unlike the true tensor of second derivatives, its approximation is
- * not necessarily symmetric. This is due to the fact that in the
- * derivation, it is not clear whether we shall consider as projected
- * second derivative the term $\nabla^2 u y_{KK'}$ or $y_{KK'}^T
- * \nabla^2 u$. Depending on which choice we take, we obtain one
- * approximation of the tensor of second derivatives or its
- * transpose. To avoid this ambiguity, as result we take the
- * symmetrized form, which is the mean value of the approximation and
- * its transpose.
- *
- * The returned value on each cell is the spectral norm of the
- * approximated tensor of second derivatives, i.e. the largest
- * eigenvalue by absolute value. This equals the largest curvature of
- * the finite element field at each cell, and the spectral norm is the
- * matrix norm associated to the $l_2$ vector norm.
- *
- * Even higher than the second derivative can be obtained along the
- * same lines as exposed above.
- *
- *
- * @sect2{Refinement indicators based on the derivatives}
- *
- * If you would like to base a refinement criterion upon these
- * approximation of the derivatives, you will have to scale the results
- * of this class by an appropriate power of the mesh width. For
- * example, since
- * $\|u-u_h\|^2_{L_2} \le C h^2 \|\nabla u\|^2_{L_2}$, it might be the
- * right thing to scale the indicators as $\eta_K = h \|\nabla u\|_K$,
- * i.e. $\eta_K = h^{1+d/2} \|\nabla u\|_{\infty;K}$, i.e. the right
- * power is $1+d/2$.
- *
- * Likewise, for the second derivative, one should choose a power of
- * the mesh size $h$ one higher than for the gradient.
- *
- *
- * @sect2{Implementation}
- *
- * The formulae for the computation of approximations to the gradient
- * and to the tensor of second derivatives shown above are very much
- * alike. The basic difference is that in one case the finite
- * difference quotiont is a scalar, while in the other case it is a
- * vector. For higher derivatives, this would be a tensor of even
- * higher rank. We then have to form the outer product of this
- * difference quotient with the distance vector $y_{KK'}$, symmetrize
- * it, contract it with the matrix $Y^{-1}$ and compute its norm. To
- * make the implementation simpler and to allow for code reuse, all
- * these operations that are dependent on the actual order of the
- * derivatives to be approximated, as well as the computation of the
- * quantities entering the difference quotient, have been separated
- * into auxiliary nested classes (names @p{Gradient} and
- * @p{SecondDerivative}) and the main algorithm is simply passed one
- * or the other data types and asks them to perform the order
- * dependent operations. The main framework that is independent of
- * this, such as finding all active neighbors, or setting up the
- * matrix $Y$ is done in the main function @p{approximate}.
- *
- * Due to this way of operation, the class may be easily extended for
- * higher oder derivatives than are presently implemented. Basically,
- * only an additional class along the lines of the derivative
- * descriptor classes @p{Gradient} and @p{SecondDerivative} has to be
- * implemented, with the respective typedefs and functions replaced by
- * the appropriate analogues for the derivative that is to be
- * approximated.
- *
- * @author Wolfgang Bangerth, 2000
- */
-class DerivativeApproximation
-{
- public:
- /**
- * This function is used to
- * obtain an approximation of the
- * gradient. Pass it the DoF
- * handler object that describes
- * the finite element field, a
- * nodal value vector, and
- * receive the cell-wise
- * Euclidian norm of the
- * approximated gradient.
- */
- template <int dim>
- static void
- approximate_gradient (const DoFHandler<dim> &dof,
- const Vector<double> &solution,
- Vector<float> &derivative_norm);
-
- /**
- * This function is the analogue
- * to the one above, computing
- * finite difference
- * approximations of the tensor
- * of second derivatives. Pass it
- * the DoF handler object that
- * describes the finite element
- * field, a nodal value vector,
- * and receive the cell-wise
- * spectral norm of the
- * approximated tensor of second
- * derivatives. The spectral norm
- * is the matrix norm associated
- * to the $l_2$ vector norm.
- */
- template <int dim>
- static void
- approximate_second_derivative (const DoFHandler<dim> &dof,
- const Vector<double> &solution,
- Vector<float> &derivative_norm);
-
- /**
- * Exception
- */
- DeclException2 (ExcInvalidVectorLength,
- int, int,
- << "Vector has length " << arg1 << ", but should have "
- << arg2);
- /**
- * Exception
- */
- DeclException0 (ExcInsufficientDirections);
-
- private:
-
- /**
- * The following class is used to
- * describe the data needed to
- * compute the finite difference
- * approximation to the gradient
- * on a cell. See the general
- * documentation of this class
- * for more information on
- * implementational details.
- *
- * @author Wolfgang Bangerth, 2000
- */
- template <int dim>
- class Gradient
- {
- public:
- /**
- * Declare which data fields have
- * to be updated for the function
- * @p{get_projected_derivative}
- * to work.
- */
- static const UpdateFlags update_flags = update_values;
-
- /**
- * Declare the data type which
- * holds the derivative described
- * by this class.
- */
- typedef Tensor<1,dim> Derivative;
-
- /**
- * Likewise declare the data type
- * that holds the derivative
- * projected to a certain
- * directions.
- */
- typedef double ProjectedDerivative;
-
- /**
- * Given an @p{FEValues} object
- * initialized to a cell, and a
- * solution vector, extract the
- * desired derivative at the
- * first quadrature point (which
- * is the only one, as we only
- * evaluate the finite element
- * field at the center of each
- * cell).
- */
- static ProjectedDerivative
- get_projected_derivative (const FEValues<dim> &fe_values,
- const Vector<double> &solution);
-
- /**
- * Return the norm of the
- * derivative object. Here, for
- * the gradient, we choose the
- * Euclidian norm of the gradient
- * vector.
- */
- static double derivative_norm (const Derivative &d);
-
- /**
- * If for the present derivative
- * order, symmetrization of the
- * derivative tensor is
- * necessary, then do so on the
- * argument.
- *
- * For the first derivatives, no
- * such thing is necessary, so
- * this function is a no-op.
- */
- static void symmetrize (Derivative &derivative_tensor);
- };
-
-
-
- /**
- * The following class is used to
- * describe the data needed to
- * compute the finite difference
- * approximation to the second
- * derivatives on a cell. See the
- * general documentation of this
- * class for more information on
- * implementational details.
- *
- * @author Wolfgang Bangerth, 2000
- */
- template <int dim>
- class SecondDerivative
- {
- public:
- /**
- * Declare which data fields have
- * to be updated for the function
- * @p{get_projected_derivative}
- * to work.
- */
- static const UpdateFlags update_flags = update_gradients;
-
- /**
- * Declare the data type which
- * holds the derivative described
- * by this class.
- */
- typedef Tensor<2,dim> Derivative;
-
- /**
- * Likewise declare the data type
- * that holds the derivative
- * projected to a certain
- * directions.
- */
- typedef Tensor<1,dim> ProjectedDerivative;
-
- /**
- * Given an @p{FEValues} object
- * initialized to a cell, and a
- * solution vector, extract the
- * desired derivative at the
- * first quadrature point (which
- * is the only one, as we only
- * evaluate the finite element
- * field at the center of each
- * cell).
- */
- static ProjectedDerivative
- get_projected_derivative (const FEValues<dim> &fe_values,
- const Vector<double> &solution);
-
- /**
- * Return the norm of the
- * derivative object. Here, for
- * the (symmetric) tensor of
- * second derivatives, we choose
- * the absolute value of the
- * largest eigenvalue, which is
- * the matrix norm associated to
- * the $l_2$ norm of vectors. It
- * is also the largest value of
- * the curvature of the solution.
- */
- static double derivative_norm (const Derivative &d);
-
- /**
- * If for the present derivative
- * order, symmetrization of the
- * derivative tensor is
- * necessary, then do so on the
- * argument.
- *
- * For the second derivatives,
- * each entry of the tensor is
- * set to the mean of its value
- * and the value of the transpose
- * element.
- *
- * Note that this function
- * actually modifies its
- * argument.
- */
- static void symmetrize (Derivative &derivative_tensor);
- };
-
- /**
- * Convenience typedef denoting
- * the range of indices on which
- * a certain thread shall
- * operate.
- */
- typedef pair<unsigned int,unsigned int> IndexInterval;
-
- /**
- * Kind of the main function of
- * this class. It is called by
- * the public entry points to
- * this class with the correct
- * template first argument and
- * then simply calls the
- * @p{approximate} function,
- * after setting up several
- * threads and doing some
- * administration that is
- * independent of the actual
- * derivative to be computed.
- */
- template <class DerivativeDescription, int dim>
- static void
- approximate_derivative (const DoFHandler<dim> &dof,
- const Vector<double> &solution,
- Vector<float> &derivative_norm);
-
- /**
- * Compute the derivative
- * approximation on the cells in
- * the range given by the third
- * parameter.
- */
- template <class DerivativeDescription, int dim>
- static void
- approximate (const DoFHandler<dim> &dof,
- const Vector<double> &solution,
- const IndexInterval &index_interval,
- Vector<float> &derivative_norm);
-};
-
-
-#endif
-
-
+++ /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 <base/thread_management.h>
-#include <base/multithread_info.h>
-#include <lac/vector.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_handler.h>
-#include <fe/fe.h>
-#include <fe/fe_values.h>
-#include <numerics/gradient_estimator.h>
-
-
-template <typename T>
-static T sqr (const T t)
-{
- return t*t;
-};
-
-
-
-
-
-template <int dim>
-inline
-typename DerivativeApproximation::Gradient<dim>::ProjectedDerivative
-DerivativeApproximation::Gradient<dim>::
-get_projected_derivative (const FEValues<dim> &fe_values,
- const Vector<double> &solution)
-{
- vector<ProjectedDerivative> values (1);
- fe_values.get_function_values (solution, values);
- return values[0];
-};
-
-
-
-template <int dim>
-inline
-double
-DerivativeApproximation::Gradient<dim>::derivative_norm (const Derivative &d)
-{
- double s = 0;
- for (unsigned int i=0; i<dim; ++i)
- s += d[i]*d[i];
- return sqrt(s);
-};
-
-
-
-template <int dim>
-inline
-void
-DerivativeApproximation::Gradient<dim>::symmetrize (Derivative &)
-{
- // nothing to do here
-};
-
-
-
-template <int dim>
-inline
-typename DerivativeApproximation::SecondDerivative<dim>::ProjectedDerivative
-DerivativeApproximation::SecondDerivative<dim>::
-get_projected_derivative (const FEValues<dim> &fe_values,
- const Vector<double> &solution)
-{
- vector<ProjectedDerivative> values (1);
- fe_values.get_function_grads (solution, values);
- return values[0];
-};
-
-
-
-template <>
-inline
-double
-DerivativeApproximation::SecondDerivative<1>::
-derivative_norm (const Derivative &d)
-{
- return fabs (d[0][0]);
-};
-
-
-
-template <>
-inline
-double
-DerivativeApproximation::SecondDerivative<2>::
-derivative_norm (const Derivative &d)
-{
- // note that d should be a
- // symmetric 2x2 tensor, so the
- // eigenvalues are:
- //
- // 1/2(a+b\pm\sqrt((a-b)^2+4c^2))
- //
- // if the d_11=a, d_22=b,
- // d_12=d_21=c
- const double radicand = sqr(d[0][0] - d[1][1]) + 4*sqr(d[0][1]);
- const double eigenvalues[2]
- = { 0.5*(d[0][0] + d[1][1] + sqrt(radicand)),
- 0.5*(d[0][0] + d[1][1] - sqrt(radicand)) };
-
- return max (fabs (eigenvalues[0]),
- fabs (eigenvalues[1]));
-};
-
-
-
-template <int dim>
-inline
-double
-DerivativeApproximation::SecondDerivative<dim>::
-derivative_norm (const Derivative &d)
-{
- // computing the spectral norm is
- // not so simple in general. it is
- // feasible for dim==3, since then
- // there are still closed form
- // expressions of the roots of the
- // third order characteristic
- // polynomial, and they can easily
- // be computed using
- // maple. however, for higher
- // dimensions, some other method
- // needs to be employed.
- Assert (false, ExcNotImplemented());
- return 0;
-};
-
-
-
-template <int dim>
-inline
-void
-DerivativeApproximation::SecondDerivative<dim>::symmetrize (Derivative &d)
-{
- // symmetrize non-diagonal entries
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=i+1; j<dim; ++j)
- {
- const double s = (d[i][j] + d[j][i]) / 2;
- d[i][j] = d[j][i] = s;
- };
-};
-
-
-
-
-template <int dim>
-void
-DerivativeApproximation::
-approximate_gradient (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &derivative_norm)
-{
- approximate_derivative<Gradient<dim>,dim> (dof_handler,
- solution,
- derivative_norm);
-};
-
-
-
-template <int dim>
-void
-DerivativeApproximation::
-approximate_second_derivative (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &derivative_norm)
-{
- approximate_derivative<SecondDerivative<dim>,dim> (dof_handler,
- solution,
- derivative_norm);
-};
-
-
-
-template <class DerivativeDescription, int dim>
-void
-DerivativeApproximation::
-approximate_derivative (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &derivative_norm)
-{
- Assert (derivative_norm.size() == dof_handler.get_tria().n_active_cells(),
- ExcInvalidVectorLength (derivative_norm.size(),
- dof_handler.get_tria().n_active_cells()));
- Assert (dof_handler.get_fe().n_components() == 1,
- ExcInternalError());
-
- const unsigned int n_threads = multithread_info.n_default_threads;
- vector<IndexInterval> index_intervals
- = Threads::split_interval (0, dof_handler.get_tria().n_active_cells(),
- n_threads);
- Threads::ThreadManager thread_manager;
- for (unsigned int i=0; i<n_threads; ++i)
- Threads::spawn (thread_manager,
- Threads::encapsulate
- (&DerivativeApproximation::
- template approximate<DerivativeDescription,dim>)
- .collect_args (dof_handler, solution,
- index_intervals[i],
- derivative_norm));
- thread_manager.wait ();
-};
-
-
-
-template <class DerivativeDescription, int dim>
-void
-DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
- const IndexInterval &index_interval,
- Vector<float> &derivative_norm)
-{
- QMidpoint<dim> midpoint_rule;
- FEValues<dim> fe_midpoint_value (dof_handler.get_fe(),
- midpoint_rule,
- UpdateFlags(DerivativeDescription::update_flags |
- update_q_points));
-
- // matrix Y=sum_i y_i y_i^T
- Tensor<2,dim> Y;
-
- // iterators over all cells and the
- // respective entries in the output
- // vector:
- Vector<float>::iterator
- derivative_norm_on_this_cell
- = derivative_norm.begin() + index_interval.first;
-
- typename DoFHandler<dim>::active_cell_iterator cell, endc;
- cell = endc = dof_handler.begin_active();
- // (static_cast to avoid warnings
- // about unsigned always >=0)
- advance (cell, static_cast<int>(index_interval.first));
- advance (endc, static_cast<int>(index_interval.second));
-
- // vector to hold iterators to all
- // active neighbors of a cell
- // reserve the maximal number of
- // active neighbors
- vector<typename DoFHandler<dim>::active_cell_iterator> active_neighbors;
- active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
- GeometryInfo<dim>::subfaces_per_face);
-
- for (; cell!=endc; ++cell, ++derivative_norm_on_this_cell)
- {
- Y.clear ();
- // vector
- // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
- // or related type for higher
- // derivatives
- typename DerivativeDescription::Derivative projected_derivative;
-
- // reinit fe values object...
- fe_midpoint_value.reinit (cell);
-
- // ...and get the value of the
- // projected derivative...
- const typename DerivativeDescription::ProjectedDerivative
- this_midpoint_value
- = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
- solution);
- // ...and the place where it lives
- const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
-
-
- // loop over all neighbors and
- // accumulate the difference
- // quotients from them. note
- // that things get a bit more
- // complicated if the neighbor
- // is more refined than the
- // present one
- //
- // to make processing simpler,
- // first collect all neighbor
- // cells in a vector, and then
- // collect the data from them
- active_neighbors.clear ();
- for (unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
- if (! cell->at_boundary(n))
- {
- typename DoFHandler<dim>::cell_iterator
- neighbor = cell->neighbor(n);
- if (neighbor->active())
- active_neighbors.push_back (neighbor);
- else
- {
- // check children
- // of
- // neighbor. note
- // that in 1d
- // children of
- // the neighbor
- // may be further
- // refined, while
- // they can't in
- // more than one
- // dimension. however,
- // in 1d the case
- // is simpler
- // since we know
- // what children
- // bound to the
- // present cell
- if (dim == 1)
- {
- typename DoFHandler<dim>::cell_iterator
- neighbor_child = neighbor;
- while (neighbor_child->has_children())
- neighbor_child = neighbor_child->child (n==0 ? 1 : 0);
-
- Assert (neighbor_child->neighbor(n==0 ? 1 : 0)==cell,
- ExcInternalError());
-
- active_neighbors.push_back (neighbor_child);
- }
- else
- // this neighbor has
- // children. find out
- // which border to the
- // present cell
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (neighbor->child(c)->neighbor(f) == cell)
- active_neighbors.push_back (neighbor->child(c));
- };
- };
-
- // now loop over all active
- // neighbors and collect the
- // data we need
- typename vector<typename DoFHandler<dim>::active_cell_iterator>::const_iterator
- neighbor_ptr = active_neighbors.begin();
- for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
- {
- const typename DoFHandler<dim>::active_cell_iterator
- neighbor = *neighbor_ptr;
-
- // reinit fe values object...
- fe_midpoint_value.reinit (neighbor);
-
- // ...and get the value of the
- // solution...
- const typename DerivativeDescription::ProjectedDerivative
- neighbor_midpoint_value
- = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
- solution);
-
- // ...and the place where it lives
- const Point<dim>
- neighbor_center = fe_midpoint_value.quadrature_point(0);
-
-
- // vector for the
- // normalized
- // direction between
- // the centers of two
- // cells
- Point<dim> y = neighbor_center - this_center;
- const double distance = sqrt(y.square());
- // normalize y
- y /= distance;
- // *** note that unlike in
- // the docs, y denotes the
- // normalized vector
- // connecting the centers
- // of the two cells, rather
- // than the normal
- // difference! ***
-
- // add up the
- // contribution of
- // this cell to Y
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- Y[i][j] += y[i] * y[j];
-
- // then update the sum
- // of difference
- // quotients
- typename DerivativeDescription::ProjectedDerivative
- projected_finite_difference
- = (neighbor_midpoint_value -
- this_midpoint_value);
- projected_finite_difference /= distance;
-
- typename DerivativeDescription::Derivative projected_derivative_update;
- outer_product (projected_derivative_update,
- y,
- projected_finite_difference);
- projected_derivative += projected_derivative_update;
- };
-
- // can we determine an
- // approximation of the
- // gradient for the present
- // cell? if so, then we need to
- // have passed over vectors y_i
- // which span the whole space,
- // otherwise we would not have
- // all components of the
- // gradient
- AssertThrow (determinant(Y) != 0,
- ExcInsufficientDirections());
-
- // first symmetrize g
- DerivativeDescription::symmetrize (projected_derivative);
-
- // compute Y^-1 g
- typename DerivativeDescription::Derivative derivative;
- Tensor<2,dim> Y_inverse = invert(Y);
-
- contract (derivative, Y_inverse, projected_derivative);
-
- *derivative_norm_on_this_cell
- = DerivativeDescription::derivative_norm (derivative);
- };
-};
-
-
-
-
-// explicit instantiations
-template
-void
-DerivativeApproximation::
-approximate_gradient (const DoFHandler<deal_II_dimension> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &derivative_norm);
-
-template
-void
-DerivativeApproximation::
-approximate_second_derivative (const DoFHandler<deal_II_dimension> &dof_handler,
- const Vector<double> &solution,
- Vector<float> &derivative_norm);
-
-
-
projects / dealii.git / commitdiff