derivative_approximation.cc

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2000 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/work_stream.h>
 
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping.h>
 
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/reference_cell.h>
#include <deal.II/grid/tria_iterator.h>
 
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/q_collection.h>
 
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/petsc_block_vector.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/trilinos_epetra_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_tpetra_block_vector.h>
#include <deal.II/lac/trilinos_tpetra_vector.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
 
#include <deal.II/numerics/derivative_approximation.h>
 
#include <algorithm>
#include <cmath>
 
DEAL_II_NAMESPACE_OPEN
 
 
 
namespace
{
  template <typename T>
  T
  sqr(const T t)
  {
    return t * t;
  }
} // namespace
 
// --- First the classes and functions that describe individual derivatives ---
 
namespace DerivativeApproximation
{
  namespace internal
  {
    template <int dim>
    class Gradient
    {
    public:
      static const UpdateFlags update_flags;
 
      using Derivative = Tensor<1, dim>;
 
      using ProjectedDerivative = Tensor<0, dim>;
 
      template <class InputVector, int spacedim>
      static ProjectedDerivative
      get_projected_derivative(const FEValues<dim, spacedim> &fe_values,
                               const InputVector             &solution,
                               const unsigned int             component);
 
      static double
      derivative_norm(const Derivative &d);
 
      static void
      symmetrize(Derivative &derivative_tensor);
    };
 
    // static variables
    template <int dim>
    const UpdateFlags Gradient<dim>::update_flags = update_values;
 
 
    template <int dim>
    template <class InputVector, int spacedim>
    typename Gradient<dim>::ProjectedDerivative
    Gradient<dim>::get_projected_derivative(
      const FEValues<dim, spacedim> &fe_values,
      const InputVector             &solution,
      const unsigned int             component)
    {
      if (fe_values.get_fe().n_components() == 1)
        {
          std::vector<typename InputVector::value_type> values(1);
          fe_values.get_function_values(solution, values);
          return values[0];
        }
      else
        {
          std::vector<Vector<typename InputVector::value_type>> values(
            1,
            Vector<typename InputVector::value_type>(
              fe_values.get_fe().n_components()));
          fe_values.get_function_values(solution, values);
          return values[0](component);
        }
    }
 
 
 
    template <int dim>
    double
    Gradient<dim>::derivative_norm(const Derivative &d)
    {
      double s = 0;
      for (unsigned int i = 0; i < dim; ++i)
        s += d[i] * d[i];
      return std::sqrt(s);
    }
 
 
 
    template <int dim>
    void
    Gradient<dim>::symmetrize(Derivative &)
    {
      // nothing to do here
    }
 
 
 
    template <int dim>
    class SecondDerivative
    {
    public:
      static const UpdateFlags update_flags;
 
      using Derivative = Tensor<2, dim>;
 
      using ProjectedDerivative = Tensor<1, dim>;
 
      template <class InputVector, int spacedim>
      static ProjectedDerivative
      get_projected_derivative(const FEValues<dim, spacedim> &fe_values,
                               const InputVector             &solution,
                               const unsigned int             component);
 
      static double
      derivative_norm(const Derivative &d);
 
      static void
      symmetrize(Derivative &derivative_tensor);
    };
 
    template <int dim>
    const UpdateFlags SecondDerivative<dim>::update_flags = update_gradients;
 
 
    template <int dim>
    template <class InputVector, int spacedim>
    typename SecondDerivative<dim>::ProjectedDerivative
    SecondDerivative<dim>::get_projected_derivative(
      const FEValues<dim, spacedim> &fe_values,
      const InputVector             &solution,
      const unsigned int             component)
    {
      if (fe_values.get_fe().n_components() == 1)
        {
          std::vector<Tensor<1, dim, typename InputVector::value_type>> values(
            1);
          fe_values.get_function_gradients(solution, values);
          return ProjectedDerivative(values[0]);
        }
      else
        {
          std::vector<
            std::vector<Tensor<1, dim, typename InputVector::value_type>>>
            values(
              1,
              std::vector<Tensor<1, dim, typename InputVector::value_type>>(
                fe_values.get_fe().n_components()));
          fe_values.get_function_gradients(solution, values);
          return ProjectedDerivative(values[0][component]);
        }
    }
 
 
 
    template <>
    double
    SecondDerivative<1>::derivative_norm(const Derivative &d)
    {
      return std::fabs(d[0][0]);
    }
 
 
 
    template <>
    double
    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] + std::sqrt(radicand)),
        0.5 * (d[0][0] + d[1][1] - std::sqrt(radicand))};
 
      return std::max(std::fabs(eigenvalues[0]), std::fabs(eigenvalues[1]));
    }
 
 
 
    template <>
    double
    SecondDerivative<3>::derivative_norm(const Derivative &d)
    {
      /*
      compute the three eigenvalues of the tensor @p{d} and take the
      largest. one could use the following maple script to generate C
      code:
 
      with(linalg);
      readlib(C);
      A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]);
      E:=eigenvals(A);
      EE:=vector(3,[E[1],E[2],E[3]]);
      C(EE);
 
      Unfortunately, with both optimized and non-optimized output, at some
      places the code `sqrt(-1.0)' is emitted, and I don't know what
      Maple intends to do with it. This happens both with Maple4 and
      Maple5.
 
      Fortunately, Roger Young provided the following Fortran code, which
      is transcribed below to C. The code uses an algorithm that uses the
      invariants of a symmetric matrix. (The translated algorithm is
      augmented by a test for R>0, since R==0 indicates that all three
      eigenvalues are equal.)
 
 
          PROGRAM MAIN
 
      C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX
      C (ROGER YOUNG, 2001)
 
          IMPLICIT NONE
 
          REAL*8 A11, A12, A13, A22, A23, A33
          REAL*8 I1, J2, J3, AM
          REAL*8 S11, S12, S13, S22, S23, S33
          REAL*8 SS12, SS23, SS13
          REAL*8 R,R3, XX,YY, THETA
          REAL*8 A1,A2,A3
          REAL*8 PI
          PARAMETER (PI=3.141592653587932384D0)
          REAL*8 A,B,C, TOL
          PARAMETER (TOL=1.D-14)
 
      C DEFINE A TEST MATRIX
 
          A11 = -1.D0
          A12 = 5.D0
          A13 = 3.D0
          A22 = -2.D0
          A23 = 0.5D0
          A33 = 4.D0
 
 
          I1 = A11 + A22 + A33
          AM = I1/3.D0
 
          S11 = A11 - AM
          S22 = A22 - AM
          S33 = A33 - AM
          S12 = A12
          S13 = A13
          S23 = A23
 
          SS12 = S12*S12
          SS23 = S23*S23
          SS13 = S13*S13
 
          J2 = S11*S11 + S22*S22 + S33*S33
          J2 = J2 + 2.D0*(SS12 + SS23 + SS13)
          J2 = J2/2.D0
 
          J3 = S11**3 + S22**3 + S33**3
          J3 = J3 + 3.D0*S11*(SS12 + SS13)
          J3 = J3 + 3.D0*S22*(SS12 + SS23)
          J3 = J3 + 3.D0*S33*(SS13 + SS23)
          J3 = J3 + 6.D0*S12*S23*S13
          J3 = J3/3.D0
 
          R = SQRT(4.D0*J2/3.D0)
          R3 = R*R*R
          XX = 4.D0*J3/R3
 
          YY = 1.D0 - DABS(XX)
          IF(YY.LE.0.D0)THEN
             IF(YY.GT.(-TOL))THEN
                WRITE(6,*)'Equal roots: XX= ',XX
                A = -(XX/DABS(XX))*SQRT(J2/3.D0)
                B = AM + A
                C = AM - 2.D0*A
                WRITE(6,*)B,' (twice) ',C
                STOP
             ELSE
                WRITE(6,*)'Error: XX= ',XX
                STOP
             ENDIF
          ENDIF
 
          THETA = (ACOS(XX))/3.D0
 
          A1 = AM + R*COS(THETA)
          A2 = AM + R*COS(THETA + 2.D0*PI/3.D0)
          A3 = AM + R*COS(THETA + 4.D0*PI/3.D0)
 
          WRITE(6,*)A1,A2,A3
 
          STOP
          END
 
       */
 
      const double am = trace(d) / 3.;
 
      // s := d - trace(d) I
      Tensor<2, 3> s = d;
      for (unsigned int i = 0; i < 3; ++i)
        s[i][i] -= am;
 
      const double ss01 = s[0][1] * s[0][1], ss12 = s[1][2] * s[1][2],
                   ss02 = s[0][2] * s[0][2];
 
      const double J2 = (s[0][0] * s[0][0] + s[1][1] * s[1][1] +
                         s[2][2] * s[2][2] + 2 * (ss01 + ss02 + ss12)) /
                        2.;
      const double J3 =
        (Utilities::fixed_power<3>(s[0][0]) +
         Utilities::fixed_power<3>(s[1][1]) +
         Utilities::fixed_power<3>(s[2][2]) + 3. * s[0][0] * (ss01 + ss02) +
         3. * s[1][1] * (ss01 + ss12) + 3. * s[2][2] * (ss02 + ss12) +
         6. * s[0][1] * s[0][2] * s[1][2]) /
        3.;
 
      const double R = std::sqrt(4. * J2 / 3.);
 
      double EE[3] = {0, 0, 0};
      // the eigenvalues are away from
      // @p{am} in the order of R. thus,
      // if R<<AM, then we have the
      // degenerate case with three
      // identical eigenvalues. check
      // this first
      if (R <= 1e-14 * std::fabs(am))
        EE[0] = EE[1] = EE[2] = am;
      else
        {
          // at least two eigenvalues are
          // distinct
          const double R3 = R * R * R;
          const double XX = 4. * J3 / R3;
          const double YY = 1. - std::fabs(XX);
 
          Assert(YY > -1e-14, ExcInternalError());
 
          if (YY < 0)
            {
              // two roots are equal
              const double a = (XX > 0 ? -1. : 1.) * R / 2;
              EE[0] = EE[1] = am + a;
              EE[2]         = am - 2. * a;
            }
          else
            {
              const double theta = std::acos(XX) / 3.;
              EE[0]              = am + R * std::cos(theta);
              EE[1] = am + R * std::cos(theta + 2. / 3. * numbers::PI);
              EE[2] = am + R * std::cos(theta + 4. / 3. * numbers::PI);
            }
        }
 
      return std::max({std::fabs(EE[0]), std::fabs(EE[1]), std::fabs(EE[2])});
    }
 
 
 
    template <int dim>
    double
    SecondDerivative<dim>::derivative_norm(const Derivative &)
    {
      // computing the spectral norm is
      // not so simple in general. it is
      // feasible for dim==3 as shown
      // above, since then there are
      // still closed form expressions of
      // the roots of the characteristic
      // polynomial, and they can easily
      // be computed using
      // maple. however, for higher
      // dimensions, some other method
      // needs to be employed. maybe some
      // steps of the power method would
      // suffice?
      DEAL_II_NOT_IMPLEMENTED();
      return 0;
    }
 
 
 
    template <int dim>
    void
    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>
    class ThirdDerivative
    {
    public:
      static const UpdateFlags update_flags;
 
      using Derivative = Tensor<3, dim>;
 
      using ProjectedDerivative = Tensor<2, dim>;
 
      template <class InputVector, int spacedim>
      static ProjectedDerivative
      get_projected_derivative(const FEValues<dim, spacedim> &fe_values,
                               const InputVector             &solution,
                               const unsigned int             component);
 
      static double
      derivative_norm(const Derivative &d);
 
      static void
      symmetrize(Derivative &derivative_tensor);
    };
 
    template <int dim>
    const UpdateFlags ThirdDerivative<dim>::update_flags = update_hessians;
 
 
    template <int dim>
    template <class InputVector, int spacedim>
    typename ThirdDerivative<dim>::ProjectedDerivative
    ThirdDerivative<dim>::get_projected_derivative(
      const FEValues<dim, spacedim> &fe_values,
      const InputVector             &solution,
      const unsigned int             component)
    {
      if (fe_values.get_fe().n_components() == 1)
        {
          std::vector<Tensor<2, dim, typename InputVector::value_type>> values(
            1);
          fe_values.get_function_hessians(solution, values);
          return ProjectedDerivative(values[0]);
        }
      else
        {
          std::vector<
            std::vector<Tensor<2, dim, typename InputVector::value_type>>>
            values(
              1,
              std::vector<Tensor<2, dim, typename InputVector::value_type>>(
                fe_values.get_fe().n_components()));
          fe_values.get_function_hessians(solution, values);
          return ProjectedDerivative(values[0][component]);
        }
    }
 
 
 
    template <>
    double
    ThirdDerivative<1>::derivative_norm(const Derivative &d)
    {
      return std::fabs(d[0][0][0]);
    }
 
 
 
    template <int dim>
    double
    ThirdDerivative<dim>::derivative_norm(const Derivative &d)
    {
      // return the Frobenius-norm. this is a
      // member function of Tensor<rank_,dim>
      return d.norm();
    }
 
 
    template <int dim>
    void
    ThirdDerivative<dim>::symmetrize(Derivative &d)
    {
      // symmetrize non-diagonal entries
 
      // first do it in the case, that i,j,k are
      // pairwise different (which can onlky happen
      // in dim >= 3)
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int j = i + 1; j < dim; ++j)
          for (unsigned int k = j + 1; k < dim; ++k)
            {
              const double s = (d[i][j][k] + d[i][k][j] + d[j][i][k] +
                                d[j][k][i] + d[k][i][j] + d[k][j][i]) /
                               6;
              d[i][j][k] = d[i][k][j] = d[j][i][k] = d[j][k][i] = d[k][i][j] =
                d[k][j][i]                                      = s;
            }
      // now do the case, where two indices are
      // equal
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int j = i + 1; j < dim; ++j)
          {
            // case 1: index i (lower one) is
            // double
            const double s = (d[i][i][j] + d[i][j][i] + d[j][i][i]) / 3;
            d[i][i][j] = d[i][j][i] = d[j][i][i] = s;
 
            // case 2: index j (higher one) is
            // double
            const double t = (d[i][j][j] + d[j][i][j] + d[j][j][i]) / 3;
            d[i][j][j] = d[j][i][j] = d[j][j][i] = t;
          }
    }
 
 
    template <int order, int dim>
    class DerivativeSelector
    {
    public:
      using DerivDescr = void;
    };
 
    template <int dim>
    class DerivativeSelector<1, dim>
    {
    public:
      using DerivDescr = Gradient<dim>;
    };
 
    template <int dim>
    class DerivativeSelector<2, dim>
    {
    public:
      using DerivDescr = SecondDerivative<dim>;
    };
 
    template <int dim>
    class DerivativeSelector<3, dim>
    {
    public:
      using DerivDescr = ThirdDerivative<dim>;
    };
  } // namespace internal
} // namespace DerivativeApproximation
 
// Dummy structures and dummy function used for WorkStream
namespace DerivativeApproximation
{
  namespace internal
  {
    namespace Assembler
    {
      struct Scratch
      {
        Scratch() = default;
      };
 
      struct CopyData
      {
        CopyData() = default;
      };
    } // namespace Assembler
  }   // namespace internal
} // namespace DerivativeApproximation
 
// --------------- now for the functions that do the actual work --------------
 
namespace DerivativeApproximation
{
  namespace internal
  {
    template <class DerivativeDescription,
              int dim,
              class InputVector,
              int spacedim>
    void
    approximate_cell(
      const Mapping<dim, spacedim>    &mapping,
      const DoFHandler<dim, spacedim> &dof_handler,
      const InputVector               &solution,
      const unsigned int               component,
      const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
      typename DerivativeDescription::Derivative &derivative)
    {
      const QMidpoint<dim> midpoint_rule;
 
      // create collection objects from
      // single quadratures, mappings,
      // and finite elements. if we have
      // an hp-DoFHandler,
      // dof_handler.get_fe() returns a
      // collection of which we do a
      // shallow copy instead
      const hp::QCollection<dim>   q_collection(midpoint_rule);
      const hp::FECollection<dim> &fe_collection =
        dof_handler.get_fe_collection();
      const hp::MappingCollection<dim> mapping_collection(mapping);
 
      hp::FEValues<dim> x_fe_midpoint_value(
        mapping_collection,
        fe_collection,
        q_collection,
        DerivativeDescription::update_flags | update_quadrature_points);
 
      // matrix Y=sum_i y_i y_i^T
      Tensor<2, dim> Y;
 
 
      // vector to hold iterators to all
      // active neighbors of a cell
      // reserve the maximal number of
      // active neighbors
      std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
        active_neighbors;
 
      active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
                               GeometryInfo<dim>::max_children_per_face);
 
      // 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...
      x_fe_midpoint_value.reinit(cell);
      const FEValues<dim> &fe_midpoint_value =
        x_fe_midpoint_value.get_present_fe_values();
 
      // ...and get the value of the
      // projected derivative...
      const typename DerivativeDescription::ProjectedDerivative
        this_midpoint_value =
          DerivativeDescription::get_projected_derivative(fe_midpoint_value,
                                                          solution,
                                                          component);
      // ...and the place where it lives
      // This needs to be a copy. If it was a reference, it would be changed
      // after the next `reinit` call of the FEValues object. clang-tidy
      // complains about this not being a reference, so we suppress the warning.
      const Point<dim> this_center =
        fe_midpoint_value.quadrature_point(0); // NOLINT
 
      // 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
      GridTools::get_active_neighbors<DoFHandler<dim, spacedim>>(
        cell, active_neighbors);
 
      // now loop over all active
      // neighbors and collect the
      // data we need
      auto neighbor_ptr = active_neighbors.begin();
      for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
        {
          const auto &neighbor = *neighbor_ptr;
 
          // reinit FE values object...
          x_fe_midpoint_value.reinit(neighbor);
          const FEValues<dim> &neighbor_fe_midpoint_value =
            x_fe_midpoint_value.get_present_fe_values();
 
          // ...and get the value of the
          // solution...
          const typename DerivativeDescription::ProjectedDerivative
            neighbor_midpoint_value =
              DerivativeDescription::get_projected_derivative(
                neighbor_fe_midpoint_value, solution, component);
 
          // ...and the place where it lives
          const Point<dim> &neighbor_center =
            neighbor_fe_midpoint_value.quadrature_point(0);
 
 
          // vector for the
          // normalized
          // direction between
          // the centers of two
          // cells
 
          Tensor<1, dim> y        = neighbor_center - this_center;
          const double   distance = y.norm();
          // 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;
 
          projected_derivative += outer_product(y, projected_finite_difference);
        }
 
      // 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());
 
      // compute Y^-1 g
      const Tensor<2, dim> Y_inverse = invert(Y);
 
      derivative = Y_inverse * projected_derivative;
 
      // finally symmetrize the derivative
      DerivativeDescription::symmetrize(derivative);
    }
 
 
 
    template <class DerivativeDescription,
              int dim,
              class InputVector,
              int spacedim>
    void
    approximate(
      const SynchronousIterators<
        std::tuple<typename DoFHandler<dim, spacedim>::active_cell_iterator,
                   Vector<float>::iterator>> &cell,
      const Mapping<dim, spacedim>           &mapping,
      const DoFHandler<dim, spacedim>        &dof_handler,
      const InputVector                      &solution,
      const unsigned int                      component)
    {
      // if the cell is not locally owned, then there is nothing to do
      if (std::get<0>(*cell)->is_locally_owned() == false)
        *std::get<1>(*cell) = 0;
      else
        {
          typename DerivativeDescription::Derivative derivative;
          // call the function doing the actual
          // work on this cell
          approximate_cell<DerivativeDescription, dim, InputVector, spacedim>(
            mapping,
            dof_handler,
            solution,
            component,
            std::get<0>(*cell),
            derivative);
 
          // evaluate the norm and fill the vector
          //*derivative_norm_on_this_cell
          *std::get<1>(*cell) =
            DerivativeDescription::derivative_norm(derivative);
        }
    }
 
 
    template <class DerivativeDescription,
              int dim,
              class InputVector,
              int spacedim>
    void
    approximate_derivative(const Mapping<dim, spacedim>    &mapping,
                           const DoFHandler<dim, spacedim> &dof_handler,
                           const InputVector               &solution,
                           const unsigned int               component,
                           Vector<float>                   &derivative_norm)
    {
      Assert(derivative_norm.size() ==
               dof_handler.get_triangulation().n_active_cells(),
             ExcVectorLengthVsNActiveCells(
               derivative_norm.size(),
               dof_handler.get_triangulation().n_active_cells()));
      AssertIndexRange(component, dof_handler.get_fe(0).n_components());
 
      using Iterators =
        std::tuple<typename DoFHandler<dim, spacedim>::active_cell_iterator,
                   Vector<float>::iterator>;
      SynchronousIterators<Iterators> begin(
        Iterators(dof_handler.begin_active(), derivative_norm.begin())),
        end(Iterators(dof_handler.end(), derivative_norm.end()));
 
      // Use type erasure to prevent workstream from seeing this function's type
      // (which in turn prevents it from instantiating for every possible
      // InputVector type)
      std::function<void(const SynchronousIterators<Iterators> &cell,
                         const Assembler::Scratch &,
                         Assembler::CopyData &)>
        worker = [&mapping,
                  &dof_handler,
                  &solution,
                  component](const SynchronousIterators<Iterators> &cell,
                             const Assembler::Scratch &,
                             Assembler::CopyData &) {
          approximate<DerivativeDescription, dim, InputVector, spacedim>(
            cell, mapping, dof_handler, solution, component);
        };
 
      // There is no need for a copier because there is no conflict between
      // threads to write in derivative_norm. Scratch and CopyData are also
      // useless.
      WorkStream::run(
        begin,
        end,
        worker,
        std::function<void(const internal::Assembler::CopyData &)>(),
        internal::Assembler::Scratch(),
        internal::Assembler::CopyData());
    }
 
  } // namespace internal
 
} // namespace DerivativeApproximation
 
 
// ------------------------ finally for the public interface of this namespace
 
namespace DerivativeApproximation
{
  template <int dim, class InputVector, int spacedim>
  void
  approximate_gradient(const Mapping<dim, spacedim>    &mapping,
                       const DoFHandler<dim, spacedim> &dof_handler,
                       const InputVector               &solution,
                       Vector<float>                   &derivative_norm,
                       const unsigned int               component)
  {
    internal::approximate_derivative<internal::Gradient<dim>, dim>(
      mapping, dof_handler, solution, component, derivative_norm);
  }
 
 
  template <int dim, class InputVector, int spacedim>
  void
  approximate_gradient(const DoFHandler<dim, spacedim> &dof_handler,
                       const InputVector               &solution,
                       Vector<float>                   &derivative_norm,
                       const unsigned int               component)
  {
    Assert(!dof_handler.get_triangulation().is_mixed_mesh(),
           ExcNotImplemented());
    const auto reference_cell =
      dof_handler.get_triangulation().get_reference_cells()[0];
    internal::approximate_derivative<internal::Gradient<dim>, dim>(
      reference_cell.template get_default_linear_mapping<spacedim>(),
      dof_handler,
      solution,
      component,
      derivative_norm);
  }
 
 
  template <int dim, class InputVector, int spacedim>
  void
  approximate_second_derivative(const Mapping<dim, spacedim>    &mapping,
                                const DoFHandler<dim, spacedim> &dof_handler,
                                const InputVector               &solution,
                                Vector<float>     &derivative_norm,
                                const unsigned int component)
  {
    internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
      mapping, dof_handler, solution, component, derivative_norm);
  }
 
 
  template <int dim, class InputVector, int spacedim>
  void
  approximate_second_derivative(const DoFHandler<dim, spacedim> &dof_handler,
                                const InputVector               &solution,
                                Vector<float>     &derivative_norm,
                                const unsigned int component)
  {
    Assert(!dof_handler.get_triangulation().is_mixed_mesh(),
           ExcNotImplemented());
    const auto reference_cell =
      dof_handler.get_triangulation().get_reference_cells()[0];
    internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
      reference_cell.template get_default_linear_mapping<spacedim>(),
      dof_handler,
      solution,
      component,
      derivative_norm);
  }
 
 
  template <int dim, int spacedim, class InputVector, int order>
  void
  approximate_derivative_tensor(
    const Mapping<dim, spacedim>    &mapping,
    const DoFHandler<dim, spacedim> &dof,
    const InputVector               &solution,
#ifndef _MSC_VER
    const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
#else
    const TriaActiveIterator<::DoFCellAccessor<dim, spacedim, false>>
      &cell,
#endif
    Tensor<order, dim> &derivative,
    const unsigned int  component)
  {
    internal::approximate_cell<
      typename internal::DerivativeSelector<order, dim>::DerivDescr>(
      mapping, dof, solution, component, cell, derivative);
  }
 
 
 
  template <int dim, int spacedim, class InputVector, int order>
  void
  approximate_derivative_tensor(
    const DoFHandler<dim, spacedim> &dof,
    const InputVector               &solution,
#ifndef _MSC_VER
    const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
#else
    const TriaActiveIterator<::DoFCellAccessor<dim, spacedim, false>>
      &cell,
#endif
    Tensor<order, dim> &derivative,
    const unsigned int  component)
  {
    // just call the respective function with Q1 mapping
    approximate_derivative_tensor(
      cell->reference_cell().template get_default_linear_mapping<spacedim>(),
      dof,
      solution,
      cell,
      derivative,
      component);
  }
 
 
 
  template <int dim, int order>
  double
  derivative_norm(const Tensor<order, dim> &derivative)
  {
    return internal::DerivativeSelector<order, dim>::DerivDescr::
      derivative_norm(derivative);
  }
 
} // namespace DerivativeApproximation
 
 
// --------------------------- explicit instantiations ---------------------
#include "numerics/derivative_approximation.inst"
 
 
DEAL_II_NAMESPACE_CLOSE