fe_values.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_fe_values_h
#define dealii_fe_values_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/point.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/vector_slice.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/table.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/mapping.h>

#include <algorithm>
#include <memory>
#include <type_traits>


// dummy include in order to have the
// definition of PetscScalar available
// without including other PETSc stuff
#ifdef DEAL_II_WITH_PETSC
#  include <petsc.h>
#endif

DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim=dim> class FEValuesBase;

namespace internal
{
  template <int dim, class NumberType=double>
  struct CurlType;

  template <class NumberType>
  struct CurlType<1, NumberType>
  {
    typedef Tensor<1,1,NumberType>     type;
  };

  template <class NumberType>
  struct CurlType<2,NumberType>
  {
    typedef Tensor<1,1,NumberType>     type;
  };

  template <class NumberType>
  struct CurlType<3,NumberType>
  {
    typedef Tensor<1,3,NumberType>     type;
  };
}




namespace FEValuesViews
{
  template <int dim, int spacedim=dim>
  class Scalar
  {
  public:
    typedef double        value_type;

    typedef ::Tensor<1,spacedim> gradient_type;

    typedef ::Tensor<2,spacedim> hessian_type;

    typedef ::Tensor<3,spacedim> third_derivative_type;

    template <typename Number>
    struct OutputType
    {
      typedef typename ProductType<Number, typename Scalar<dim,spacedim>::value_type>::type value_type;

      typedef typename ProductType<Number, typename Scalar<dim,spacedim>::gradient_type>::type gradient_type;

      typedef typename ProductType<Number, typename Scalar<dim,spacedim>::value_type>::type laplacian_type;

      typedef typename ProductType<Number, typename Scalar<dim,spacedim>::hessian_type>::type hessian_type;

      typedef typename ProductType<Number, typename Scalar<dim,spacedim>::third_derivative_type>::type third_derivative_type;
    };

    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component;

      unsigned int row_index;
    };

    Scalar ();

    Scalar (const FEValuesBase<dim,spacedim> &fe_values_base,
            const unsigned int       component);

    Scalar &operator= (const Scalar<dim,spacedim> &);

    value_type
    value (const unsigned int shape_function,
           const unsigned int q_point) const;

    gradient_type
    gradient (const unsigned int shape_function,
              const unsigned int q_point) const;

    hessian_type
    hessian (const unsigned int shape_function,
             const unsigned int q_point) const;

    third_derivative_type
    third_derivative (const unsigned int shape_function,
                      const unsigned int q_point) const;

    template <class InputVector>
    void get_function_values (const InputVector &fe_function,
                              std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const;

    template <class InputVector>
    void get_function_values_from_local_dof_values (const InputVector &dof_values,
                                                    std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const;

    template <class InputVector>
    void get_function_gradients (const InputVector &fe_function,
                                 std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const;

    template <class InputVector>
    void get_function_gradients_from_local_dof_values (const InputVector &dof_values,
                                                       std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const;

    template <class InputVector>
    void get_function_hessians (const InputVector &fe_function,
                                std::vector<typename ProductType<hessian_type,typename InputVector::value_type>::type> &hessians) const;

    template <class InputVector>
    void get_function_hessians_from_local_dof_values (const InputVector &dof_values,
                                                      std::vector<typename OutputType<typename InputVector::value_type>::hessian_type> &hessians) const;


    template <class InputVector>
    void get_function_laplacians (const InputVector &fe_function,
                                  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &laplacians) const;

    template <class InputVector>
    void get_function_laplacians_from_local_dof_values (const InputVector &dof_values,
                                                        std::vector<typename OutputType<typename InputVector::value_type>::laplacian_type> &laplacians) const;


    template <class InputVector>
    void get_function_third_derivatives (const InputVector &fe_function,
                                         std::vector<typename ProductType<third_derivative_type,
                                         typename InputVector::value_type>::type> &third_derivatives) const;

    template <class InputVector>
    void get_function_third_derivatives_from_local_dof_values (const InputVector &dof_values,
                                                               std::vector<typename OutputType<typename InputVector::value_type>::third_derivative_type> &third_derivatives) const;


  private:
    const SmartPointer<const FEValuesBase<dim,spacedim> > fe_values;

    const unsigned int component;

    std::vector<ShapeFunctionData> shape_function_data;
  };



  template <int dim, int spacedim=dim>
  class Vector
  {
  public:
    typedef ::Tensor<1,spacedim>          value_type;

    typedef ::Tensor<2,spacedim>          gradient_type;

    typedef ::SymmetricTensor<2,spacedim> symmetric_gradient_type;

    typedef double                 divergence_type;

    typedef typename ::internal::CurlType<spacedim>::type   curl_type;

    typedef ::Tensor<3,spacedim>          hessian_type;

    typedef ::Tensor<4,spacedim>          third_derivative_type;

    template <typename Number>
    struct OutputType
    {
      typedef typename ProductType<Number, typename Vector<dim,spacedim>::value_type>::type value_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::gradient_type>::type gradient_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::symmetric_gradient_type>::type symmetric_gradient_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::divergence_type>::type divergence_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::value_type>::type laplacian_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::curl_type>::type curl_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::hessian_type>::type hessian_type;

      typedef typename ProductType<Number, typename Vector<dim,spacedim>::third_derivative_type>::type third_derivative_type;
    };

    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component[spacedim];

      unsigned int row_index[spacedim];

      int          single_nonzero_component;
      unsigned int single_nonzero_component_index;
    };

    Vector ();

    Vector (const FEValuesBase<dim,spacedim> &fe_values_base,
            const unsigned int first_vector_component);

    Vector &operator= (const Vector<dim,spacedim> &);

    value_type
    value (const unsigned int shape_function,
           const unsigned int q_point) const;

    gradient_type
    gradient (const unsigned int shape_function,
              const unsigned int q_point) const;

    symmetric_gradient_type
    symmetric_gradient (const unsigned int shape_function,
                        const unsigned int q_point) const;

    divergence_type
    divergence (const unsigned int shape_function,
                const unsigned int q_point) const;

    curl_type
    curl (const unsigned int shape_function,
          const unsigned int q_point) const;

    hessian_type
    hessian (const unsigned int shape_function,
             const unsigned int q_point) const;

    third_derivative_type
    third_derivative (const unsigned int shape_function,
                      const unsigned int q_point) const;

    template <class InputVector>
    void get_function_values (const InputVector &fe_function,
                              std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const;

    template <class InputVector>
    void get_function_values_from_local_dof_values (const InputVector &dof_values,
                                                    std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const;

    template <class InputVector>
    void get_function_gradients (const InputVector &fe_function,
                                 std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const;

    template <class InputVector>
    void get_function_gradients_from_local_dof_values (const InputVector &dof_values,
                                                       std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const;

    template <class InputVector>
    void
    get_function_symmetric_gradients (const InputVector &fe_function,
                                      std::vector<typename ProductType<symmetric_gradient_type,typename InputVector::value_type>::type> &symmetric_gradients) const;

    template <class InputVector>
    void
    get_function_symmetric_gradients_from_local_dof_values (const InputVector &dof_values,
                                                            std::vector<typename OutputType<typename InputVector::value_type>::symmetric_gradient_type> &symmetric_gradients) const;

    template <class InputVector>
    void get_function_divergences (const InputVector &fe_function,
                                   std::vector<typename ProductType<divergence_type,typename InputVector::value_type>::type> &divergences) const;

    template <class InputVector>
    void get_function_divergences_from_local_dof_values (const InputVector &dof_values,
                                                         std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const;

    template <class InputVector>
    void get_function_curls (const InputVector &fe_function,
                             std::vector<typename ProductType<curl_type,typename InputVector::value_type>::type> &curls) const;

    template <class InputVector>
    void get_function_curls_from_local_dof_values (const InputVector &dof_values,
                                                   std::vector<typename OutputType<typename InputVector::value_type>::curl_type> &curls) const;

    template <class InputVector>
    void get_function_hessians (const InputVector &fe_function,
                                std::vector<typename ProductType<hessian_type,typename InputVector::value_type>::type> &hessians) const;

    template <class InputVector>
    void get_function_hessians_from_local_dof_values (const InputVector &dof_values,
                                                      std::vector<typename OutputType<typename InputVector::value_type>::hessian_type> &hessians) const;

    template <class InputVector>
    void get_function_laplacians (const InputVector &fe_function,
                                  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &laplacians) const;

    template <class InputVector>
    void get_function_laplacians_from_local_dof_values (const InputVector &dof_values,
                                                        std::vector<typename OutputType<typename InputVector::value_type>::laplacian_type> &laplacians) const;

    template <class InputVector>
    void get_function_third_derivatives (const InputVector &fe_function,
                                         std::vector<typename ProductType<third_derivative_type,
                                         typename InputVector::value_type>::type> &third_derivatives) const;

    template <class InputVector>
    void get_function_third_derivatives_from_local_dof_values (const InputVector &dof_values,
                                                               std::vector<typename OutputType<typename InputVector::value_type>::third_derivative_type> &third_derivatives) const;

  private:
    const SmartPointer<const FEValuesBase<dim,spacedim> > fe_values;

    const unsigned int first_vector_component;

    std::vector<ShapeFunctionData> shape_function_data;
  };


  template <int rank, int dim, int spacedim = dim>
  class SymmetricTensor;

  template <int dim, int spacedim>
  class SymmetricTensor<2,dim,spacedim>
  {
  public:
    typedef ::SymmetricTensor<2, spacedim> value_type;

    typedef ::Tensor<1, spacedim> divergence_type;

    template <typename Number>
    struct OutputType
    {
      typedef typename ProductType<Number, typename SymmetricTensor<2,dim,spacedim>::value_type>::type value_type;

      typedef typename ProductType<Number, typename SymmetricTensor<2,dim,spacedim>::divergence_type>::type divergence_type;
    };

    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component[value_type::n_independent_components];

      unsigned int row_index[value_type::n_independent_components];

      int single_nonzero_component;

      unsigned int single_nonzero_component_index;
    };

    SymmetricTensor();

    SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values_base,
                    const unsigned int first_tensor_component);

    SymmetricTensor &operator=(const SymmetricTensor<2, dim, spacedim> &);

    value_type
    value (const unsigned int shape_function,
           const unsigned int q_point) const;

    divergence_type
    divergence (const unsigned int shape_function,
                const unsigned int q_point) const;

    template <class InputVector>
    void get_function_values (const InputVector &fe_function,
                              std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const;

    template <class InputVector>
    void get_function_values_from_local_dof_values (const InputVector &dof_values,
                                                    std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const;

    template <class InputVector>
    void get_function_divergences (const InputVector &fe_function,
                                   std::vector<typename ProductType<divergence_type,typename InputVector::value_type>::type> &divergences) const;

    template <class InputVector>
    void get_function_divergences_from_local_dof_values (const InputVector &dof_values,
                                                         std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const;

  private:
    const SmartPointer<const FEValuesBase<dim, spacedim> > fe_values;

    const unsigned int first_tensor_component;

    std::vector<ShapeFunctionData> shape_function_data;
  };


  template <int rank, int dim, int spacedim = dim>
  class Tensor;

  template <int dim, int spacedim>
  class Tensor<2,dim,spacedim>
  {
  public:

    typedef ::Tensor<2, spacedim> value_type;

    typedef ::Tensor<1, spacedim> divergence_type;

    typedef ::Tensor<3, spacedim> gradient_type;

    template <typename Number>
    struct OutputType
    {
      typedef typename ProductType<Number, typename Tensor<2,dim,spacedim>::value_type>::type value_type;

      typedef typename ProductType<Number, typename Tensor<2,dim,spacedim>::divergence_type>::type divergence_type;

      typedef typename ProductType<Number, typename Tensor<2,dim,spacedim>::gradient_type>::type gradient_type;
    };

    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component[value_type::n_independent_components];

      unsigned int row_index[value_type::n_independent_components];

      int single_nonzero_component;

      unsigned int single_nonzero_component_index;
    };

    Tensor();

    Tensor(const FEValuesBase<dim, spacedim> &fe_values_base,
           const unsigned int first_tensor_component);


    Tensor &operator=(const Tensor<2, dim, spacedim> &);

    value_type
    value (const unsigned int shape_function,
           const unsigned int q_point) const;

    divergence_type
    divergence (const unsigned int shape_function,
                const unsigned int q_point) const;

    gradient_type
    gradient (const unsigned int shape_function,
              const unsigned int q_point) const;

    template <class InputVector>
    void get_function_values (const InputVector &fe_function,
                              std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const;

    template <class InputVector>
    void get_function_values_from_local_dof_values (const InputVector &dof_values,
                                                    std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const;

    template <class InputVector>
    void get_function_divergences (const InputVector &fe_function,
                                   std::vector<typename ProductType<divergence_type,typename InputVector::value_type>::type> &divergences) const;

    template <class InputVector>
    void get_function_divergences_from_local_dof_values (const InputVector &dof_values,
                                                         std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const;

    template <class InputVector>
    void get_function_gradients (const InputVector &fe_function,
                                 std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const;

    template <class InputVector>
    void get_function_gradients_from_local_dof_values (const InputVector &dof_values,
                                                       std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const;

  private:
    const SmartPointer<const FEValuesBase<dim, spacedim> > fe_values;

    const unsigned int first_tensor_component;

    std::vector<ShapeFunctionData> shape_function_data;
  };

}


namespace internal
{
  namespace FEValuesViews
  {
    template <int dim, int spacedim>
    struct Cache
    {
      std::vector<::FEValuesViews::Scalar<dim,spacedim> > scalars;
      std::vector<::FEValuesViews::Vector<dim,spacedim> > vectors;
      std::vector<::FEValuesViews::SymmetricTensor<2,dim,spacedim> >
      symmetric_second_order_tensors;
      std::vector<::FEValuesViews::Tensor<2,dim,spacedim> >
      second_order_tensors;

      Cache (const FEValuesBase<dim,spacedim> &fe_values);
    };
  }
}



template <int dim, int spacedim>
class FEValuesBase :
  public Subscriptor
{
public:
  static const unsigned int dimension = dim;

  static const unsigned int space_dimension = spacedim;

  const unsigned int n_quadrature_points;

  const unsigned int dofs_per_cell;


  FEValuesBase (const unsigned int n_q_points,
                const unsigned int dofs_per_cell,
                const UpdateFlags update_flags,
                const Mapping<dim,spacedim> &mapping,
                const FiniteElement<dim,spacedim> &fe);


  ~FEValuesBase ();




  const double &shape_value (const unsigned int function_no,
                             const unsigned int point_no) const;

  double shape_value_component (const unsigned int function_no,
                                const unsigned int point_no,
                                const unsigned int component) const;

  const Tensor<1,spacedim> &
  shape_grad (const unsigned int function_no,
              const unsigned int quadrature_point) const;

  Tensor<1,spacedim>
  shape_grad_component (const unsigned int function_no,
                        const unsigned int point_no,
                        const unsigned int component) const;

  const Tensor<2,spacedim> &
  shape_hessian (const unsigned int function_no,
                 const unsigned int point_no) const;

  Tensor<2,spacedim>
  shape_hessian_component (const unsigned int function_no,
                           const unsigned int point_no,
                           const unsigned int component) const;

  const Tensor<3,spacedim> &
  shape_3rd_derivative (const unsigned int function_no,
                        const unsigned int point_no) const;

  Tensor<3,spacedim>
  shape_3rd_derivative_component (const unsigned int function_no,
                                  const unsigned int point_no,
                                  const unsigned int component) const;



  template <class InputVector>
  void get_function_values (const InputVector &fe_function,
                            std::vector<typename InputVector::value_type> &values) const;

  template <class InputVector>
  void get_function_values (const InputVector       &fe_function,
                            std::vector<Vector<typename InputVector::value_type> > &values) const;

  template <class InputVector>
  void get_function_values (const InputVector &fe_function,
                            const VectorSlice<const std::vector<types::global_dof_index> > &indices,
                            std::vector<typename InputVector::value_type> &values) const;

  template <class InputVector>
  void get_function_values (const InputVector &fe_function,
                            const VectorSlice<const std::vector<types::global_dof_index> > &indices,
                            std::vector<Vector<typename InputVector::value_type> > &values) const;


  template <class InputVector>
  void get_function_values (const InputVector &fe_function,
                            const VectorSlice<const std::vector<types::global_dof_index> > &indices,
                            VectorSlice<std::vector<std::vector<typename InputVector::value_type> > > values,
                            const bool quadrature_points_fastest) const;



  template <class InputVector>
  void get_function_gradients (const InputVector      &fe_function,
                               std::vector<Tensor<1,spacedim,typename InputVector::value_type> > &gradients) const;

  template <class InputVector>
  void get_function_gradients (const InputVector               &fe_function,
                               std::vector<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > > &gradients) const;

  template <class InputVector>
  void get_function_gradients (const InputVector &fe_function,
                               const VectorSlice<const std::vector<types::global_dof_index> > &indices,
                               std::vector<Tensor<1,spacedim,typename InputVector::value_type> > &gradients) const;

  template <class InputVector>
  void get_function_gradients (const InputVector &fe_function,
                               const VectorSlice<const std::vector<types::global_dof_index> > &indices,
                               VectorSlice<std::vector<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > > > gradients,
                               bool quadrature_points_fastest = false) const;



  template <class InputVector>
  void
  get_function_hessians (const InputVector &fe_function,
                         std::vector<Tensor<2,spacedim,typename InputVector::value_type> > &hessians) const;

  template <class InputVector>
  void
  get_function_hessians (const InputVector      &fe_function,
                         std::vector<std::vector<Tensor<2,spacedim,typename InputVector::value_type> > > &hessians,
                         bool quadrature_points_fastest = false) const;

  template <class InputVector>
  void get_function_hessians (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<Tensor<2,spacedim,typename InputVector::value_type> > &hessians) const;

  template <class InputVector>
  void get_function_hessians (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    VectorSlice<std::vector<std::vector<Tensor<2,spacedim,typename InputVector::value_type> > > > hessians,
    bool quadrature_points_fastest = false) const;

  template <class InputVector>
  void
  get_function_laplacians (const InputVector &fe_function,
                           std::vector<typename InputVector::value_type> &laplacians) const;

  template <class InputVector>
  void
  get_function_laplacians (const InputVector      &fe_function,
                           std::vector<Vector<typename InputVector::value_type> > &laplacians) const;

  template <class InputVector>
  void get_function_laplacians (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<typename InputVector::value_type> &laplacians) const;

  template <class InputVector>
  void get_function_laplacians (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<Vector<typename InputVector::value_type> > &laplacians) const;

  template <class InputVector>
  void get_function_laplacians (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<std::vector<typename InputVector::value_type> > &laplacians,
    bool quadrature_points_fastest = false) const;



  template <class InputVector>
  void
  get_function_third_derivatives (const InputVector &fe_function,
                                  std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const;

  template <class InputVector>
  void
  get_function_third_derivatives (const InputVector      &fe_function,
                                  std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > &third_derivatives,
                                  bool quadrature_points_fastest = false) const;

  template <class InputVector>
  void get_function_third_derivatives (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const;

  template <class InputVector>
  void get_function_third_derivatives (
    const InputVector &fe_function,
    const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    VectorSlice<std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > > third_derivatives,
    bool quadrature_points_fastest = false) const;



  const Point<spacedim> &
  quadrature_point (const unsigned int q) const;

  const std::vector<Point<spacedim> > &get_quadrature_points () const;

  double JxW (const unsigned int quadrature_point) const;

  const std::vector<double> &get_JxW_values () const;

  const DerivativeForm<1,dim,spacedim> &jacobian (const unsigned int quadrature_point) const;

  const std::vector<DerivativeForm<1,dim,spacedim> > &get_jacobians () const;

  const DerivativeForm<2,dim,spacedim> &jacobian_grad (const unsigned int quadrature_point) const;

  const std::vector<DerivativeForm<2,dim,spacedim> > &get_jacobian_grads () const;

  const Tensor<3,spacedim> &jacobian_pushed_forward_grad (const unsigned int quadrature_point) const;

  const std::vector<Tensor<3,spacedim> > &get_jacobian_pushed_forward_grads () const;

  const DerivativeForm<3,dim,spacedim> &jacobian_2nd_derivative (const unsigned int quadrature_point) const;

  const std::vector<DerivativeForm<3,dim,spacedim> > &get_jacobian_2nd_derivatives () const;

  const Tensor<4,spacedim> &jacobian_pushed_forward_2nd_derivative (const unsigned int quadrature_point) const;

  const std::vector<Tensor<4,spacedim> > &get_jacobian_pushed_forward_2nd_derivatives () const;

  const DerivativeForm<4,dim,spacedim> &jacobian_3rd_derivative (const unsigned int quadrature_point) const;

  const std::vector<DerivativeForm<4,dim,spacedim> > &get_jacobian_3rd_derivatives () const;

  const Tensor<5,spacedim> &jacobian_pushed_forward_3rd_derivative (const unsigned int quadrature_point) const;

  const std::vector<Tensor<5,spacedim> > &get_jacobian_pushed_forward_3rd_derivatives () const;

  const DerivativeForm<1,spacedim,dim> &inverse_jacobian (const unsigned int quadrature_point) const;

  const std::vector<DerivativeForm<1,spacedim,dim> > &get_inverse_jacobians () const;

  const Tensor<1,spacedim> &normal_vector (const unsigned int i) const;

  DEAL_II_DEPRECATED
  const std::vector<Tensor<1,spacedim> > &get_all_normal_vectors () const;

  const std::vector<Tensor<1,spacedim> > &get_normal_vectors () const;




  const FEValuesViews::Scalar<dim,spacedim> &
  operator[] (const FEValuesExtractors::Scalar &scalar) const;

  const FEValuesViews::Vector<dim,spacedim> &
  operator[] (const FEValuesExtractors::Vector &vector) const;

  const FEValuesViews::SymmetricTensor<2,dim,spacedim> &
  operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const;


  const FEValuesViews::Tensor<2,dim,spacedim> &
  operator[] (const FEValuesExtractors::Tensor<2> &tensor) const;




  const Mapping<dim,spacedim> &get_mapping () const;

  const FiniteElement<dim,spacedim> &get_fe () const;

  UpdateFlags get_update_flags () const;

  const typename Triangulation<dim,spacedim>::cell_iterator get_cell () const;

  CellSimilarity::Similarity get_cell_similarity () const;

  std::size_t memory_consumption () const;


  DeclException1 (ExcAccessToUninitializedField,
                  char *,
                  << "You are requesting information from an FEValues/FEFaceValues/FESubfaceValues "
                  << "object for which this kind of information has not been computed. What "
                  << "information these objects compute is determined by the update_* flags you "
                  << "pass to the constructor. Here, the operation you are attempting requires "
                  << "the <" << arg1 << "> flag to be set, but it was apparently not specified "
                  << "upon construction.");

  DeclExceptionMsg (ExcFEDontMatch,
                    "The FiniteElement you provided to FEValues and the FiniteElement that belongs "
                    "to the DoFHandler that provided the cell iterator do not match.");
  DeclException1 (ExcShapeFunctionNotPrimitive,
                  int,
                  << "The shape function with index " << arg1
                  << " is not primitive, i.e. it is vector-valued and "
                  << "has more than one non-zero vector component. This "
                  << "function cannot be called for these shape functions. "
                  << "Maybe you want to use the same function with the "
                  << "_component suffix?");

  DeclExceptionMsg (ExcFENotPrimitive,
                    "The given FiniteElement is not a primitive element but the requested operation "
                    "only works for those. See FiniteElement::is_primitive() for more information.");

protected:
  class CellIteratorBase;

  template <typename CI> class CellIterator;
  class TriaCellIterator;

  std::unique_ptr<const CellIteratorBase> present_cell;

  boost::signals2::connection tria_listener_refinement;

  boost::signals2::connection tria_listener_mesh_transform;

  void invalidate_present_cell ();

  void
  maybe_invalidate_previous_present_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell);

  const SmartPointer<const Mapping<dim,spacedim>,FEValuesBase<dim,spacedim> > mapping;

  std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase> mapping_data;

  ::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> mapping_output;


  const SmartPointer<const FiniteElement<dim,spacedim>,FEValuesBase<dim,spacedim> > fe;

  std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase> fe_data;

  ::internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim> finite_element_output;


  UpdateFlags          update_flags;

  UpdateFlags compute_update_flags (const UpdateFlags update_flags) const;

  CellSimilarity::Similarity cell_similarity;

  void
  check_cell_similarity (const typename Triangulation<dim,spacedim>::cell_iterator &cell);

private:
  FEValuesBase (const FEValuesBase &);

  FEValuesBase &operator= (const FEValuesBase &);

  ::internal::FEValuesViews::Cache<dim,spacedim> fe_values_views_cache;

  template <int, int> friend class FEValuesViews::Scalar;
  template <int, int> friend class FEValuesViews::Vector;
  template <int, int, int> friend class FEValuesViews::SymmetricTensor;
  template <int, int, int> friend class FEValuesViews::Tensor;
};



template <int dim, int spacedim=dim>
class FEValues : public FEValuesBase<dim,spacedim>
{
public:
  static const unsigned int integral_dimension = dim;

  FEValues (const Mapping<dim,spacedim>       &mapping,
            const FiniteElement<dim,spacedim> &fe,
            const Quadrature<dim>             &quadrature,
            const UpdateFlags                  update_flags);

  FEValues (const FiniteElement<dim,spacedim> &fe,
            const Quadrature<dim>             &quadrature,
            const UpdateFlags                  update_flags);

  template <template <int, int> class DoFHandlerType, bool level_dof_access>
  void reinit (const TriaIterator<DoFCellAccessor<DoFHandlerType<dim,spacedim>,level_dof_access> > &cell);

  void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell);

  const Quadrature<dim> &get_quadrature () const;

  std::size_t memory_consumption () const;

  const FEValues<dim,spacedim> &get_present_fe_values () const;

private:
  const Quadrature<dim> quadrature;

  void initialize (const UpdateFlags update_flags);

  void do_reinit ();
};


template <int dim, int spacedim=dim>
class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
{
public:
  static const unsigned int integral_dimension = dim-1;

  FEFaceValuesBase (const unsigned int                 n_q_points,
                    const unsigned int                 dofs_per_cell,
                    const UpdateFlags                  update_flags,
                    const Mapping<dim,spacedim>       &mapping,
                    const FiniteElement<dim,spacedim> &fe,
                    const Quadrature<dim-1>&           quadrature);

  const Tensor<1,spacedim> &boundary_form (const unsigned int i) const;

  const std::vector<Tensor<1,spacedim> > &get_boundary_forms () const;

  unsigned int get_face_index() const;

  const Quadrature<dim-1> & get_quadrature () const;

  std::size_t memory_consumption () const;

protected:

  unsigned int present_face_index;

  const Quadrature<dim-1> quadrature;
};



template <int dim, int spacedim=dim>
class FEFaceValues : public FEFaceValuesBase<dim,spacedim>
{
public:
  static const unsigned int dimension = dim;

  static const unsigned int space_dimension = spacedim;

  static const unsigned int integral_dimension = dim-1;

  FEFaceValues (const Mapping<dim,spacedim>       &mapping,
                const FiniteElement<dim,spacedim> &fe,
                const Quadrature<dim-1>           &quadrature,
                const UpdateFlags                  update_flags);

  FEFaceValues (const FiniteElement<dim,spacedim> &fe,
                const Quadrature<dim-1>           &quadrature,
                const UpdateFlags                  update_flags);

  template <template <int, int> class DoFHandlerType, bool level_dof_access>
  void reinit (const TriaIterator<DoFCellAccessor<DoFHandlerType<dim,spacedim>,level_dof_access> > &cell,
               const unsigned int face_no);

  void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
               const unsigned int                                         face_no);

  const FEFaceValues<dim,spacedim> &get_present_fe_values () const;
private:

  void initialize (const UpdateFlags update_flags);

  void do_reinit (const unsigned int face_no);
};


template <int dim, int spacedim=dim>
class FESubfaceValues : public FEFaceValuesBase<dim,spacedim>
{
public:
  static const unsigned int dimension = dim;

  static const unsigned int space_dimension = spacedim;

  static const unsigned int integral_dimension = dim-1;

  FESubfaceValues (const Mapping<dim,spacedim>       &mapping,
                   const FiniteElement<dim,spacedim> &fe,
                   const Quadrature<dim-1>  &face_quadrature,
                   const UpdateFlags         update_flags);

  FESubfaceValues (const FiniteElement<dim,spacedim> &fe,
                   const Quadrature<dim-1>  &face_quadrature,
                   const UpdateFlags         update_flags);

  template <template <int, int> class DoFHandlerType, bool level_dof_access>
  void reinit (const TriaIterator<DoFCellAccessor<DoFHandlerType<dim,spacedim>,level_dof_access> > &cell,
               const unsigned int face_no,
               const unsigned int subface_no);

  void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
               const unsigned int                    face_no,
               const unsigned int                    subface_no);

  const FESubfaceValues<dim,spacedim> &get_present_fe_values () const;

  DeclException0 (ExcReinitCalledWithBoundaryFace);

  DeclException0 (ExcFaceHasNoSubfaces);

private:

  void initialize (const UpdateFlags update_flags);

  void do_reinit (const unsigned int face_no,
                  const unsigned int subface_no);
};


#ifndef DOXYGEN


/*------------------------ Inline functions: namespace FEValuesViews --------*/

namespace FEValuesViews
{
  template <int dim, int spacedim>
  inline
  typename Scalar<dim,spacedim>::value_type
  Scalar<dim,spacedim>::value (const unsigned int shape_function,
                               const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_values,
            ((typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values"))));

    // an adaptation of the FEValuesBase::shape_value_component function
    // except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output.shape_values(shape_function_data[shape_function]
                                                           .row_index,
                                                           q_point);
    else
      return 0;
  }



  template <int dim, int spacedim>
  inline
  typename Scalar<dim,spacedim>::gradient_type
  Scalar<dim,spacedim>::gradient (const unsigned int shape_function,
                                  const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    // an adaptation of the FEValuesBase::shape_grad_component
    // function except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function]
                                                              .row_index][q_point];
    else
      return gradient_type();
  }



  template <int dim, int spacedim>
  inline
  typename Scalar<dim,spacedim>::hessian_type
  Scalar<dim,spacedim>::hessian (const unsigned int shape_function,
                                 const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_hessians,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));

    // an adaptation of the FEValuesBase::shape_hessian_component
    // function except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output.shape_hessians[shape_function_data[shape_function].row_index][q_point];
    else
      return hessian_type();
  }



  template <int dim, int spacedim>
  inline
  typename Scalar<dim,spacedim>::third_derivative_type
  Scalar<dim,spacedim>::third_derivative (const unsigned int shape_function,
                                          const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_3rd_derivatives,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_3rd_derivatives")));

    // an adaptation of the FEValuesBase::shape_3rdderivative_component
    // function except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output.shape_3rd_derivatives[shape_function_data[shape_function].row_index][q_point];
    else
      return third_derivative_type();
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::value_type
  Vector<dim,spacedim>::value (const unsigned int shape_function,
                               const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_values,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));

    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return value_type();
    else if (snc != -1)
      {
        value_type return_value;
        return_value[shape_function_data[shape_function].single_nonzero_component_index]
          = fe_values->finite_element_output.shape_values(snc,q_point);
        return return_value;
      }
    else
      {
        value_type return_value;
        for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values->finite_element_output.shape_values(shape_function_data[shape_function].row_index[d],q_point);

        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::gradient_type
  Vector<dim,spacedim>::gradient (const unsigned int shape_function,
                                  const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return gradient_type();
    else if (snc != -1)
      {
        gradient_type return_value;
        return_value[shape_function_data[shape_function].single_nonzero_component_index]
          = fe_values->finite_element_output.shape_gradients[snc][q_point];
        return return_value;
      }
    else
      {
        gradient_type return_value;
        for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];

        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::divergence_type
  Vector<dim,spacedim>::divergence (const unsigned int shape_function,
                                    const unsigned int q_point) const
  {
    // this function works like in the case above
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return divergence_type();
    else if (snc != -1)
      return
        fe_values->finite_element_output.shape_gradients[snc][q_point][shape_function_data[shape_function].single_nonzero_component_index];
    else
      {
        divergence_type return_value = 0;
        for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value
            += fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];

        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::curl_type
  Vector<dim,spacedim>::curl (const unsigned int shape_function, const unsigned int q_point) const
  {
    // this function works like in the case above

    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;

    if (snc == -2)
      return curl_type ();

    else
      switch (dim)
        {
        case 1:
        {
          Assert (false, ExcMessage("Computing the curl in 1d is not a useful operation"));
          return curl_type ();
        }

        case 2:
        {
          if (snc != -1)
            {
              curl_type return_value;

              // the single nonzero component can only be zero or one in 2d
              if (shape_function_data[shape_function].single_nonzero_component_index == 0)
                return_value[0] = -1.0 * fe_values->finite_element_output.shape_gradients[snc][q_point][1];
              else
                return_value[0] = fe_values->finite_element_output.shape_gradients[snc][q_point][0];

              return return_value;
            }

          else
            {
              curl_type return_value;

              return_value[0] = 0.0;

              if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
                return_value[0]
                -= fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];

              if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
                return_value[0]
                += fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];

              return return_value;
            }
        }

        case 3:
        {
          if (snc != -1)
            {
              curl_type return_value;

              switch (shape_function_data[shape_function].single_nonzero_component_index)
                {
                case 0:
                {
                  return_value[0] = 0;
                  return_value[1] = fe_values->finite_element_output.shape_gradients[snc][q_point][2];
                  return_value[2] = -1.0 * fe_values->finite_element_output.shape_gradients[snc][q_point][1];
                  return return_value;
                }

                case 1:
                {
                  return_value[0] = -1.0 * fe_values->finite_element_output.shape_gradients[snc][q_point][2];
                  return_value[1] = 0;
                  return_value[2] = fe_values->finite_element_output.shape_gradients[snc][q_point][0];
                  return return_value;
                }

                default:
                {
                  return_value[0] = fe_values->finite_element_output.shape_gradients[snc][q_point][1];
                  return_value[1] = -1.0 * fe_values->finite_element_output.shape_gradients[snc][q_point][0];
                  return_value[2] = 0;
                  return return_value;
                }
                }
            }

          else
            {
              curl_type return_value;

              for (unsigned int i = 0; i < dim; ++i)
                return_value[i] = 0.0;

              if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
                {
                  return_value[1]
                  += fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][2];
                  return_value[2]
                  -= fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
                }

              if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
                {
                  return_value[0]
                  -= fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][2];
                  return_value[2]
                  += fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
                }

              if (shape_function_data[shape_function].is_nonzero_shape_function_component[2])
                {
                  return_value[0]
                  += fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][1];
                  return_value[1]
                  -= fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][0];
                }

              return return_value;
            }
        }
        }
    // should not end up here
    Assert (false, ExcInternalError());
    return curl_type();
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::hessian_type
  Vector<dim,spacedim>::hessian (const unsigned int shape_function,
                                 const unsigned int q_point) const
  {
    // this function works like in the case above
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_hessians,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));

    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return hessian_type();
    else if (snc != -1)
      {
        hessian_type return_value;
        return_value[shape_function_data[shape_function].single_nonzero_component_index]
          = fe_values->finite_element_output.shape_hessians[snc][q_point];
        return return_value;
      }
    else
      {
        hessian_type return_value;
        for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values->finite_element_output.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point];

        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::third_derivative_type
  Vector<dim,spacedim>::third_derivative (const unsigned int shape_function,
                                          const unsigned int q_point) const
  {
    // this function works like in the case above
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_3rd_derivatives,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_3rd_derivatives")));

    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return third_derivative_type();
    else if (snc != -1)
      {
        third_derivative_type return_value;
        return_value[shape_function_data[shape_function].single_nonzero_component_index]
          = fe_values->finite_element_output.shape_3rd_derivatives[snc][q_point];
        return return_value;
      }
    else
      {
        third_derivative_type return_value;
        for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values->finite_element_output.shape_3rd_derivatives[shape_function_data[shape_function].row_index[d]][q_point];

        return return_value;
      }
  }



  namespace
  {
    inline
    ::SymmetricTensor<2,1>
    symmetrize_single_row (const unsigned int n,
                           const ::Tensor<1,1> &t)
    {
      Assert (n < 1, ExcIndexRange (n, 0, 1));
      (void)n;

      const double array[1] = { t[0] };
      return ::SymmetricTensor<2,1>(array);
    }



    inline
    ::SymmetricTensor<2,2>
    symmetrize_single_row (const unsigned int n,
                           const ::Tensor<1,2> &t)
    {
      switch (n)
        {
        case 0:
        {
          const double array[3] = { t[0], 0, t[1]/2 };
          return ::SymmetricTensor<2,2>(array);
        }
        case 1:
        {
          const double array[3] = { 0, t[1], t[0]/2 };
          return ::SymmetricTensor<2,2>(array);
        }
        default:
        {
          Assert (false, ExcIndexRange (n, 0, 2));
          return ::SymmetricTensor<2,2>();
        }
        }
    }



    inline
    ::SymmetricTensor<2,3>
    symmetrize_single_row (const unsigned int n,
                           const ::Tensor<1,3> &t)
    {
      switch (n)
        {
        case 0:
        {
          const double array[6] = { t[0], 0, 0, t[1]/2, t[2]/2, 0 };
          return ::SymmetricTensor<2,3>(array);
        }
        case 1:
        {
          const double array[6] = { 0, t[1], 0, t[0]/2, 0, t[2]/2 };
          return ::SymmetricTensor<2,3>(array);
        }
        case 2:
        {
          const double array[6] = { 0, 0, t[2], 0, t[0]/2, t[1]/2 };
          return ::SymmetricTensor<2,3>(array);
        }
        default:
        {
          Assert (false, ExcIndexRange (n, 0, 3));
          return ::SymmetricTensor<2,3>();
        }
        }
    }
  }



  template <int dim, int spacedim>
  inline
  typename Vector<dim,spacedim>::symmetric_gradient_type
  Vector<dim,spacedim>::symmetric_gradient (const unsigned int shape_function,
                                            const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    // same as for the scalar case except that we have one more index
    const int snc = shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return symmetric_gradient_type();
    else if (snc != -1)
      return symmetrize_single_row (shape_function_data[shape_function].single_nonzero_component_index,
                                    fe_values->finite_element_output.shape_gradients[snc][q_point]);
    else
      {
        gradient_type return_value;
        for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values->finite_element_output.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];

        return symmetrize(return_value);
      }
  }



  template <int dim, int spacedim>
  inline
  typename SymmetricTensor<2, dim, spacedim>::value_type
  SymmetricTensor<2, dim, spacedim>::value (const unsigned int shape_function,
                                            const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_values,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));

    // similar to the vector case where we have more then one index and we need
    // to convert between unrolled and component indexing for tensors
    const int snc
      = shape_function_data[shape_function].single_nonzero_component;

    if (snc == -2)
      {
        // shape function is zero for the selected components
        return value_type();

      }
    else if (snc != -1)
      {
        value_type return_value;
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        return_value[value_type::unrolled_to_component_indices(comp)]
          = fe_values->finite_element_output.shape_values(snc,q_point);
        return return_value;
      }
    else
      {
        value_type return_value;
        for (unsigned int d = 0; d < value_type::n_independent_components; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[value_type::unrolled_to_component_indices(d)]
              = fe_values->finite_element_output.shape_values(shape_function_data[shape_function].row_index[d],q_point);
        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename SymmetricTensor<2, dim, spacedim>::divergence_type
  SymmetricTensor<2, dim, spacedim>::divergence(const unsigned int shape_function,
                                                const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    const int snc = shape_function_data[shape_function].single_nonzero_component;

    if (snc == -2)
      {
        // shape function is zero for the selected components
        return divergence_type();
      }
    else if (snc != -1)
      {
        // we have a single non-zero component when the symmetric tensor is
        // represented in unrolled form. this implies we potentially have
        // two non-zero components when represented in component form!  we
        // will only have one non-zero entry if the non-zero component lies on
        // the diagonal of the tensor.
        //
        // the divergence of a second-order tensor is a first order tensor.
        //
        // assume the second-order tensor is A with components A_{ij}.  then
        // A_{ij} = A_{ji} and there is only one (if diagonal) or two non-zero
        // entries in the tensorial representation.  define the
        // divergence as:
        // b_i := \dfrac{\partial phi_{ij}}{\partial x_j}.
        // (which is incidentally also
        // b_j := \dfrac{\partial phi_{ij}}{\partial x_i}).
        // In both cases, a sum is implied.
        //
        // Now, we know the nonzero component in unrolled form: it is indicated
        // by 'snc'. we can figure out which tensor components belong to this:
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const unsigned int ii = value_type::unrolled_to_component_indices(comp)[0];
        const unsigned int jj = value_type::unrolled_to_component_indices(comp)[1];

        // given the form of the divergence above, if ii=jj there is only a
        // single nonzero component of the full tensor and the gradient
        // equals
        // b_ii := \dfrac{\partial phi_{ii,ii}}{\partial x_ii}.
        // all other entries of 'b' are zero
        //
        // on the other hand, if ii!=jj, then there are two nonzero entries in
        // the full tensor and
        // b_ii := \dfrac{\partial phi_{ii,jj}}{\partial x_ii}.
        // b_jj := \dfrac{\partial phi_{ii,jj}}{\partial x_jj}.
        // again, all other entries of 'b' are zero
        const ::Tensor<1, spacedim> &phi_grad = fe_values->finite_element_output.shape_gradients[snc][q_point];

        divergence_type return_value;
        return_value[ii] = phi_grad[jj];

        if (ii != jj)
          return_value[jj] = phi_grad[ii];

        return return_value;

      }
    else
      {
        Assert (false, ExcNotImplemented());
        divergence_type return_value;
        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Tensor<2, dim, spacedim>::value_type
  Tensor<2, dim, spacedim>::value (const unsigned int shape_function,
                                   const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_values,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));

    // similar to the vector case where we have more then one index and we need
    // to convert between unrolled and component indexing for tensors
    const int snc
      = shape_function_data[shape_function].single_nonzero_component;

    if (snc == -2)
      {
        // shape function is zero for the selected components
        return value_type();
      }
    else if (snc != -1)
      {
        value_type return_value;
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const TableIndices<2> indices = ::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
        return_value[indices] = fe_values->finite_element_output.shape_values(snc,q_point);
        return return_value;
      }
    else
      {
        value_type return_value;
        for (unsigned int d = 0; d < dim*dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            {
              const TableIndices<2> indices = ::Tensor<2,spacedim>::unrolled_to_component_indices(d);
              return_value[indices]
                = fe_values->finite_element_output.shape_values(shape_function_data[shape_function].row_index[d],q_point);
            }
        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Tensor<2, dim, spacedim>::divergence_type
  Tensor<2, dim, spacedim>::divergence(const unsigned int shape_function,
                                       const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    const int snc = shape_function_data[shape_function].single_nonzero_component;

    if (snc == -2)
      {
        // shape function is zero for the selected components
        return divergence_type();
      }
    else if (snc != -1)
      {
        // we have a single non-zero component when the tensor is
        // represented in unrolled form.
        //
        // the divergence of a second-order tensor is a first order tensor.
        //
        // assume the second-order tensor is A with components A_{ij},
        // then divergence is d_i := \frac{\partial A_{ij}}{\partial x_j}
        //
        // Now, we know the nonzero component in unrolled form: it is indicated
        // by 'snc'. we can figure out which tensor components belong to this:
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const TableIndices<2> indices = ::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
        const unsigned int ii = indices[0];
        const unsigned int jj = indices[1];

        const ::Tensor<1, spacedim> &phi_grad = fe_values->finite_element_output.shape_gradients[snc][q_point];

        divergence_type return_value;
        // note that we contract \nabla from the right
        return_value[ii] = phi_grad[jj];

        return return_value;
      }
    else
      {
        Assert (false, ExcNotImplemented());
        divergence_type return_value;
        return return_value;
      }
  }



  template <int dim, int spacedim>
  inline
  typename Tensor<2, dim, spacedim>::gradient_type
  Tensor<2, dim, spacedim>::gradient(const unsigned int shape_function,
                                     const unsigned int q_point) const
  {
    Assert (shape_function < fe_values->fe->dofs_per_cell,
            ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
    Assert (fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));

    const int snc = shape_function_data[shape_function].single_nonzero_component;

    if (snc == -2)
      {
        // shape function is zero for the selected components
        return gradient_type();
      }
    else if (snc != -1)
      {
        // we have a single non-zero component when the tensor is
        // represented in unrolled form.
        //
        // the gradient of a second-order tensor is a third order tensor.
        //
        // assume the second-order tensor is A with components A_{ij},
        // then gradient is B_{ijk} := \frac{\partial A_{ij}}{\partial x_k}
        //
        // Now, we know the nonzero component in unrolled form: it is indicated
        // by 'snc'. we can figure out which tensor components belong to this:
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const TableIndices<2> indices = ::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
        const unsigned int ii = indices[0];
        const unsigned int jj = indices[1];

        const ::Tensor<1, spacedim> &phi_grad = fe_values->finite_element_output.shape_gradients[snc][q_point];

        gradient_type return_value;
        return_value[ii][jj] = phi_grad;

        return return_value;
      }
    else
      {
        Assert (false, ExcNotImplemented());
        gradient_type return_value;
        return return_value;
      }
  }

}



/*------------------------ Inline functions: FEValuesBase ------------------------*/



template <int dim, int spacedim>
inline
const FEValuesViews::Scalar<dim,spacedim> &
FEValuesBase<dim,spacedim>::
operator[] (const FEValuesExtractors::Scalar &scalar) const
{
  Assert (scalar.component < fe_values_views_cache.scalars.size(),
          ExcIndexRange (scalar.component,
                         0, fe_values_views_cache.scalars.size()));

  return fe_values_views_cache.scalars[scalar.component];
}



template <int dim, int spacedim>
inline
const FEValuesViews::Vector<dim,spacedim> &
FEValuesBase<dim,spacedim>::
operator[] (const FEValuesExtractors::Vector &vector) const
{
  Assert (vector.first_vector_component <
          fe_values_views_cache.vectors.size(),
          ExcIndexRange (vector.first_vector_component,
                         0, fe_values_views_cache.vectors.size()));

  return fe_values_views_cache.vectors[vector.first_vector_component];
}



template <int dim, int spacedim>
inline
const FEValuesViews::SymmetricTensor<2,dim,spacedim> &
FEValuesBase<dim,spacedim>::
operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const
{
  Assert (tensor.first_tensor_component <
          fe_values_views_cache.symmetric_second_order_tensors.size(),
          ExcIndexRange (tensor.first_tensor_component,
                         0, fe_values_views_cache.symmetric_second_order_tensors.size()));

  return fe_values_views_cache.symmetric_second_order_tensors[tensor.first_tensor_component];
}



template <int dim, int spacedim>
inline
const FEValuesViews::Tensor<2,dim,spacedim> &
FEValuesBase<dim,spacedim>::
operator[] (const FEValuesExtractors::Tensor<2> &tensor) const
{
  Assert (tensor.first_tensor_component <
          fe_values_views_cache.second_order_tensors.size(),
          ExcIndexRange (tensor.first_tensor_component,
                         0, fe_values_views_cache.second_order_tensors.size()));

  return fe_values_views_cache.second_order_tensors[tensor.first_tensor_component];
}



template <int dim, int spacedim>
inline
const double &
FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
                                         const unsigned int j) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_values,
          ExcAccessToUninitializedField("update_values"));
  Assert (fe->is_primitive (i),
          ExcShapeFunctionNotPrimitive(i));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_values(i,j);
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int
      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_values(row, j);
    }
}



template <int dim, int spacedim>
inline
double
FEValuesBase<dim,spacedim>::shape_value_component (const unsigned int i,
                                                   const unsigned int j,
                                                   const unsigned int component) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_values,
          ExcAccessToUninitializedField("update_values"));
  Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return 0;

  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int
  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_values(row, j);
}



template <int dim, int spacedim>
inline
const Tensor<1,spacedim> &
FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
                                        const unsigned int j) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_gradients,
          ExcAccessToUninitializedField("update_gradients"));
  Assert (fe->is_primitive (i),
          ExcShapeFunctionNotPrimitive(i));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_gradients[i][j];
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int
      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_gradients[row][j];
    }
}



template <int dim, int spacedim>
inline
Tensor<1,spacedim>
FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
                                                  const unsigned int j,
                                                  const unsigned int component) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_gradients,
          ExcAccessToUninitializedField("update_gradients"));
  Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return Tensor<1,spacedim>();

  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int
  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_gradients[row][j];
}



template <int dim, int spacedim>
inline
const Tensor<2,spacedim> &
FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
                                           const unsigned int j) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
  Assert (fe->is_primitive (i),
          ExcShapeFunctionNotPrimitive(i));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_hessians[i][j];
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int
      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_hessians[row][j];
    }
}



template <int dim, int spacedim>
inline
Tensor<2,spacedim>
FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
                                                     const unsigned int j,
                                                     const unsigned int component) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
  Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return Tensor<2,spacedim>();

  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int
  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_hessians[row][j];
}



template <int dim, int spacedim>
inline
const Tensor<3,spacedim> &
FEValuesBase<dim,spacedim>::shape_3rd_derivative (const unsigned int i,
                                                  const unsigned int j) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_3rd_derivatives"));
  Assert (fe->is_primitive (i),
          ExcShapeFunctionNotPrimitive(i));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_3rd_derivatives[i][j];
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int
      row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_3rd_derivatives[row][j];
    }
}



template <int dim, int spacedim>
inline
Tensor<3,spacedim>
FEValuesBase<dim,spacedim>::shape_3rd_derivative_component (const unsigned int i,
                                                            const unsigned int j,
                                                            const unsigned int component) const
{
  Assert (i < fe->dofs_per_cell,
          ExcIndexRange (i, 0, fe->dofs_per_cell));
  Assert (this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_3rd_derivatives"));
  Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return Tensor<3,spacedim>();

  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int
  row = this->finite_element_output.shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_3rd_derivatives[row][j];
}



template <int dim, int spacedim>
inline
const FiniteElement<dim,spacedim> &
FEValuesBase<dim,spacedim>::get_fe () const
{
  return *fe;
}



template <int dim, int spacedim>
inline
const Mapping<dim,spacedim> &
FEValuesBase<dim,spacedim>::get_mapping () const
{
  return *mapping;
}



template <int dim, int spacedim>
inline
UpdateFlags
FEValuesBase<dim,spacedim>::get_update_flags () const
{
  return this->update_flags;
}



template <int dim, int spacedim>
inline
const std::vector<Point<spacedim> > &
FEValuesBase<dim,spacedim>::get_quadrature_points () const
{
  Assert (this->update_flags & update_quadrature_points,
          ExcAccessToUninitializedField("update_quadrature_points"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.quadrature_points;
}



template <int dim, int spacedim>
inline
const std::vector<double> &
FEValuesBase<dim,spacedim>::get_JxW_values () const
{
  Assert (this->update_flags & update_JxW_values,
          ExcAccessToUninitializedField("update_JxW_values"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.JxW_values;
}



template <int dim, int spacedim>
inline
const std::vector<DerivativeForm<1,dim,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobians () const
{
  Assert (this->update_flags & update_jacobians,
          ExcAccessToUninitializedField("update_jacobians"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobians;
}



template <int dim, int spacedim>
inline
const std::vector<DerivativeForm<2,dim,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobian_grads () const
{
  Assert (this->update_flags & update_jacobian_grads,
          ExcAccessToUninitializedField("update_jacobians_grads"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_grads;
}



template <int dim, int spacedim>
inline
const Tensor<3,spacedim> &
FEValuesBase<dim,spacedim>::jacobian_pushed_forward_grad (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobian_pushed_forward_grads,
          ExcAccessToUninitializedField("update_jacobian_pushed_forward_grads"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_pushed_forward_grads[i];
}



template <int dim, int spacedim>
inline
const std::vector<Tensor<3,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobian_pushed_forward_grads () const
{
  Assert (this->update_flags & update_jacobian_pushed_forward_grads,
          ExcAccessToUninitializedField("update_jacobian_pushed_forward_grads"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_pushed_forward_grads;
}



template <int dim, int spacedim>
inline
const DerivativeForm<3,dim,spacedim> &
FEValuesBase<dim,spacedim>::jacobian_2nd_derivative (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobian_2nd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_2nd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_2nd_derivatives[i];
}



template <int dim, int spacedim>
inline
const std::vector<DerivativeForm<3,dim,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobian_2nd_derivatives () const
{
  Assert (this->update_flags & update_jacobian_2nd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_2nd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_2nd_derivatives;
}



template <int dim, int spacedim>
inline
const Tensor<4,spacedim> &
FEValuesBase<dim,spacedim>::jacobian_pushed_forward_2nd_derivative (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobian_pushed_forward_2nd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_pushed_forward_2nd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_pushed_forward_2nd_derivatives[i];
}



template <int dim, int spacedim>
inline
const std::vector<Tensor<4,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobian_pushed_forward_2nd_derivatives () const
{
  Assert (this->update_flags & update_jacobian_pushed_forward_2nd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_pushed_forward_2nd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_pushed_forward_2nd_derivatives;
}



template <int dim, int spacedim>
inline
const DerivativeForm<4,dim,spacedim> &
FEValuesBase<dim,spacedim>::jacobian_3rd_derivative (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobian_3rd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_3rd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_3rd_derivatives[i];
}



template <int dim, int spacedim>
inline
const std::vector<DerivativeForm<4,dim,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobian_3rd_derivatives () const
{
  Assert (this->update_flags & update_jacobian_3rd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_3rd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_3rd_derivatives;
}



template <int dim, int spacedim>
inline
const Tensor<5,spacedim> &
FEValuesBase<dim,spacedim>::jacobian_pushed_forward_3rd_derivative (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobian_pushed_forward_3rd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_pushed_forward_3rd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_pushed_forward_3rd_derivatives[i];
}



template <int dim, int spacedim>
inline
const std::vector<Tensor<5,spacedim> > &
FEValuesBase<dim,spacedim>::get_jacobian_pushed_forward_3rd_derivatives () const
{
  Assert (this->update_flags & update_jacobian_pushed_forward_3rd_derivatives,
          ExcAccessToUninitializedField("update_jacobian_pushed_forward_3rd_derivatives"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.jacobian_pushed_forward_3rd_derivatives;
}



template <int dim, int spacedim>
inline
const std::vector<DerivativeForm<1,spacedim,dim> > &
FEValuesBase<dim,spacedim>::get_inverse_jacobians () const
{
  Assert (this->update_flags & update_inverse_jacobians,
          ExcAccessToUninitializedField("update_inverse_jacobians"));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));
  return this->mapping_output.inverse_jacobians;
}



template <int dim, int spacedim>
inline
const Point<spacedim> &
FEValuesBase<dim,spacedim>::quadrature_point (const unsigned int i) const
{
  Assert (this->update_flags & update_quadrature_points,
          ExcAccessToUninitializedField("update_quadrature_points"));
  Assert (i<this->mapping_output.quadrature_points.size(),
          ExcIndexRange(i, 0, this->mapping_output.quadrature_points.size()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  return this->mapping_output.quadrature_points[i];
}



template <int dim, int spacedim>
inline
double
FEValuesBase<dim,spacedim>::JxW (const unsigned int i) const
{
  Assert (this->update_flags & update_JxW_values,
          ExcAccessToUninitializedField("update_JxW_values"));
  Assert (i<this->mapping_output.JxW_values.size(),
          ExcIndexRange(i, 0, this->mapping_output.JxW_values.size()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  return this->mapping_output.JxW_values[i];
}



template <int dim, int spacedim>
inline
const DerivativeForm<1,dim,spacedim> &
FEValuesBase<dim,spacedim>::jacobian (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobians,
          ExcAccessToUninitializedField("update_jacobians"));
  Assert (i<this->mapping_output.jacobians.size(),
          ExcIndexRange(i, 0, this->mapping_output.jacobians.size()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  return this->mapping_output.jacobians[i];
}



template <int dim, int spacedim>
inline
const DerivativeForm<2,dim,spacedim> &
FEValuesBase<dim,spacedim>::jacobian_grad (const unsigned int i) const
{
  Assert (this->update_flags & update_jacobian_grads,
          ExcAccessToUninitializedField("update_jacobians_grads"));
  Assert (i<this->mapping_output.jacobian_grads.size(),
          ExcIndexRange(i, 0, this->mapping_output.jacobian_grads.size()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  return this->mapping_output.jacobian_grads[i];
}



template <int dim, int spacedim>
inline
const DerivativeForm<1,spacedim,dim> &
FEValuesBase<dim,spacedim>::inverse_jacobian (const unsigned int i) const
{
  Assert (this->update_flags & update_inverse_jacobians,
          ExcAccessToUninitializedField("update_inverse_jacobians"));
  Assert (i<this->mapping_output.inverse_jacobians.size(),
          ExcIndexRange(i, 0, this->mapping_output.inverse_jacobians.size()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  return this->mapping_output.inverse_jacobians[i];
}



template <int dim, int spacedim>
inline
const Tensor<1,spacedim> &
FEValuesBase<dim,spacedim>::normal_vector (const unsigned int i) const
{
  Assert (this->update_flags & update_normal_vectors,
          (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_normal_vectors")));
  Assert (i<this->mapping_output.normal_vectors.size(),
          ExcIndexRange(i, 0, this->mapping_output.normal_vectors.size()));
  Assert (present_cell.get() != nullptr,
          ExcMessage ("FEValues object is not reinit'ed to any cell"));

  return this->mapping_output.normal_vectors[i];
}



/*------------------------ Inline functions: FEValues ----------------------------*/


template <int dim, int spacedim>
inline
const Quadrature<dim> &
FEValues<dim,spacedim>::get_quadrature () const
{
  return quadrature;
}



template <int dim, int spacedim>
inline
const FEValues<dim,spacedim> &
FEValues<dim,spacedim>::get_present_fe_values () const
{
  return *this;
}


/*------------------------ Inline functions: FEFaceValuesBase --------------------*/


template <int dim, int spacedim>
inline
unsigned int
FEFaceValuesBase<dim,spacedim>::get_face_index () const
{
  return present_face_index;
}


/*------------------------ Inline functions: FE*FaceValues --------------------*/

template <int dim, int spacedim>
inline
const Quadrature<dim-1> &
FEFaceValuesBase<dim,spacedim>::get_quadrature () const
{
  return quadrature;
}



template <int dim, int spacedim>
inline
const FEFaceValues<dim,spacedim> &
FEFaceValues<dim,spacedim>::get_present_fe_values () const
{
  return *this;
}



template <int dim, int spacedim>
inline
const FESubfaceValues<dim,spacedim> &
FESubfaceValues<dim,spacedim>::get_present_fe_values () const
{
  return *this;
}



template <int dim, int spacedim>
inline
const Tensor<1,spacedim> &
FEFaceValuesBase<dim,spacedim>::boundary_form (const unsigned int i) const
{
  Assert (i<this->mapping_output.boundary_forms.size(),
          ExcIndexRange(i, 0, this->mapping_output.boundary_forms.size()));
  Assert (this->update_flags & update_boundary_forms,
          (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_boundary_forms")));

  return this->mapping_output.boundary_forms[i];
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif