fe_values_views.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2023 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_fe_values_views_h
#define dealii_fe_values_views_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/lazy.h>
#include <deal.II/base/observer_pointer.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
 
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values_extractors.h>
 
#include <deal.II/lac/read_vector.h>
 
#include <type_traits>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declaration
#ifndef DOXYGEN
template <int dim, int spacedim = dim>
class FEValuesBase;
#endif
 
namespace internal
{
  template <int dim, typename NumberType = double>
  struct CurlType;
 
  template <typename NumberType>
  struct CurlType<1, NumberType>
  {
    using type = Tensor<1, 1, NumberType>;
  };
 
  template <typename NumberType>
  struct CurlType<2, NumberType>
  {
    using type = Tensor<1, 1, NumberType>;
  };
 
  template <typename NumberType>
  struct CurlType<3, NumberType>
  {
    using type = Tensor<1, 3, NumberType>;
  };
} // namespace internal
 
 
 
namespace FEValuesViews
{
  template <int dim, int spacedim = dim>
  class Scalar
  {
  public:
    using value_type = double;
 
    using gradient_type = ::Tensor<1, spacedim>;
 
    using hessian_type = ::Tensor<2, spacedim>;
 
    using third_derivative_type = ::Tensor<3, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_gradient_type =
      typename ProductType<Number, gradient_type>::type;
 
    template <typename Number>
    using solution_laplacian_type =
      typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_hessian_type =
      typename ProductType<Number, hessian_type>::type;
 
    template <typename Number>
    using solution_third_derivative_type =
      typename ProductType<Number, third_derivative_type>::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(const Scalar<dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Scalar(Scalar<dim, spacedim> &&) = default;
 
    ~Scalar() = default;
 
    Scalar &
    operator=(const Scalar<dim, spacedim> &) = delete;
 
    Scalar &
    operator=(Scalar<dim, spacedim> &&) noexcept = default;
 
    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 <typename Number>
    void
    get_function_values(const ReadVector<Number>                 &fe_function,
                        std::vector<solution_value_type<Number>> &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <typename Number>
    void
    get_function_gradients(
      const ReadVector<Number>                    &fe_function,
      std::vector<solution_gradient_type<Number>> &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <typename Number>
    void
    get_function_hessians(
      const ReadVector<Number>                   &fe_function,
      std::vector<solution_hessian_type<Number>> &hessians) const;
 
    template <class InputVector>
    void
    get_function_hessians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_hessian_type<typename InputVector::value_type>>
        &hessians) const;
 
 
    template <typename Number>
    void
    get_function_laplacians(
      const ReadVector<Number>                     &fe_function,
      std::vector<solution_laplacian_type<Number>> &laplacians) const;
 
    template <class InputVector>
    void
    get_function_laplacians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_laplacian_type<typename InputVector::value_type>>
        &laplacians) const;
 
 
    template <typename Number>
    void
    get_function_third_derivatives(
      const ReadVector<Number>                            &fe_function,
      std::vector<solution_third_derivative_type<Number>> &third_derivatives)
      const;
 
    template <class InputVector>
    void
    get_function_third_derivatives_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<
        solution_third_derivative_type<typename InputVector::value_type>>
        &third_derivatives) const;
 
 
  private:
    ObserverPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    unsigned int component;
 
    std::vector<ShapeFunctionData> shape_function_data;
  };
 
 
 
  template <int dim, int spacedim = dim>
  class Vector
  {
  public:
    using value_type = ::Tensor<1, spacedim>;
 
    using gradient_type = ::Tensor<2, spacedim>;
 
    using symmetric_gradient_type = ::SymmetricTensor<2, spacedim>;
 
    using divergence_type = double;
 
    using curl_type = typename ::internal::CurlType<spacedim>::type;
 
    using hessian_type = ::Tensor<3, spacedim>;
 
    using third_derivative_type = ::Tensor<4, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_gradient_type =
      typename ProductType<Number, gradient_type>::type;
 
    template <typename Number>
    using solution_symmetric_gradient_type =
      typename ProductType<Number, symmetric_gradient_type>::type;
 
    template <typename Number>
    using solution_divergence_type =
      typename ProductType<Number, divergence_type>::type;
 
    template <typename Number>
    using solution_laplacian_type =
      typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_curl_type = typename ProductType<Number, curl_type>::type;
 
    template <typename Number>
    using solution_hessian_type =
      typename ProductType<Number, hessian_type>::type;
 
    template <typename Number>
    using solution_third_derivative_type =
      typename ProductType<Number, third_derivative_type>::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(const Vector<dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Vector(Vector<dim, spacedim> &&) = default;
 
    ~Vector() = default;
 
    Vector &
    operator=(const Vector<dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Vector &
    operator=(Vector<dim, spacedim> &&) = default; // NOLINT
 
    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 <typename Number>
    void
    get_function_values(const ReadVector<Number>                 &fe_function,
                        std::vector<solution_value_type<Number>> &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <typename Number>
    void
    get_function_gradients(
      const ReadVector<Number>                    &fe_function,
      std::vector<solution_gradient_type<Number>> &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <typename Number>
    void
    get_function_symmetric_gradients(
      const ReadVector<Number> &fe_function,
      std::vector<solution_symmetric_gradient_type<Number>>
        &symmetric_gradients) const;
 
    template <class InputVector>
    void
    get_function_symmetric_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<
        solution_symmetric_gradient_type<typename InputVector::value_type>>
        &symmetric_gradients) const;
 
    template <typename Number>
    void
    get_function_divergences(
      const ReadVector<Number>                      &fe_function,
      std::vector<solution_divergence_type<Number>> &divergences) const;
 
    template <class InputVector>
    void
    get_function_divergences_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <typename Number>
    void
    get_function_curls(const ReadVector<Number>                &fe_function,
                       std::vector<solution_curl_type<Number>> &curls) const;
 
    template <class InputVector>
    void
    get_function_curls_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_curl_type<typename InputVector::value_type>> &curls)
      const;
 
    template <typename Number>
    void
    get_function_hessians(
      const ReadVector<Number>                   &fe_function,
      std::vector<solution_hessian_type<Number>> &hessians) const;
 
    template <class InputVector>
    void
    get_function_hessians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_hessian_type<typename InputVector::value_type>>
        &hessians) const;
 
    template <typename Number>
    void
    get_function_laplacians(
      const ReadVector<Number>                     &fe_function,
      std::vector<solution_laplacian_type<Number>> &laplacians) const;
 
    template <class InputVector>
    void
    get_function_laplacians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_laplacian_type<typename InputVector::value_type>>
        &laplacians) const;
 
    template <typename Number>
    void
    get_function_third_derivatives(
      const ReadVector<Number>                            &fe_function,
      std::vector<solution_third_derivative_type<Number>> &third_derivatives)
      const;
 
    template <class InputVector>
    void
    get_function_third_derivatives_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<
        solution_third_derivative_type<typename InputVector::value_type>>
        &third_derivatives) const;
 
  private:
    ObserverPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    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:
    using value_type = ::SymmetricTensor<2, spacedim>;
 
    using divergence_type = ::Tensor<1, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_divergence_type =
      typename ProductType<Number, divergence_type>::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(const SymmetricTensor<2, dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    SymmetricTensor(SymmetricTensor<2, dim, spacedim> &&) = default;
 
    SymmetricTensor &
    operator=(const SymmetricTensor<2, dim, spacedim> &) = delete;
 
    SymmetricTensor &
    operator=(SymmetricTensor<2, dim, spacedim> &&) noexcept = default;
 
    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 <typename Number>
    void
    get_function_values(const ReadVector<Number>                 &fe_function,
                        std::vector<solution_value_type<Number>> &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <typename Number>
    void
    get_function_divergences(
      const ReadVector<Number>                      &fe_function,
      std::vector<solution_divergence_type<Number>> &divergences) const;
 
    template <class InputVector>
    void
    get_function_divergences_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
  private:
    ObserverPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    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:
    using value_type = ::Tensor<2, spacedim>;
 
    using divergence_type = ::Tensor<1, spacedim>;
 
    using gradient_type = ::Tensor<3, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_divergence_type =
      typename ProductType<Number, divergence_type>::type;
 
    template <typename Number>
    using solution_gradient_type =
      typename ProductType<Number, gradient_type>::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 Tensor<2, dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Tensor(Tensor<2, dim, spacedim> &&) = default;
 
    ~Tensor() = default;
 
    Tensor(const FEValuesBase<dim, spacedim> &fe_values_base,
           const unsigned int                 first_tensor_component);
 
 
    Tensor &
    operator=(const Tensor<2, dim, spacedim> &) = delete;
 
    Tensor &
    operator=(Tensor<2, dim, spacedim> &&) = default; // NOLINT
 
    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 <typename Number>
    void
    get_function_values(const ReadVector<Number>                 &fe_function,
                        std::vector<solution_value_type<Number>> &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <typename Number>
    void
    get_function_divergences(
      const ReadVector<Number>                      &fe_function,
      std::vector<solution_divergence_type<Number>> &divergences) const;
 
    template <class InputVector>
    void
    get_function_divergences_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <typename Number>
    void
    get_function_gradients(
      const ReadVector<Number>                    &fe_function,
      std::vector<solution_gradient_type<Number>> &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
  private:
    ObserverPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    unsigned int first_tensor_component;
 
    std::vector<ShapeFunctionData> shape_function_data;
  };
 
} // namespace FEValuesViews
 
 
namespace internal
{
  namespace FEValuesViews
  {
    template <int dim, int spacedim, typename Extractor>
    struct ViewType
    {};
 
    template <int dim, int spacedim>
    struct ViewType<dim, spacedim, FEValuesExtractors::Scalar>
    {
      using type = typename ::FEValuesViews::Scalar<dim, spacedim>;
    };
 
    template <int dim, int spacedim>
    struct ViewType<dim, spacedim, FEValuesExtractors::Vector>
    {
      using type = typename ::FEValuesViews::Vector<dim, spacedim>;
    };
 
    template <int dim, int spacedim, int rank>
    struct ViewType<dim, spacedim, FEValuesExtractors::Tensor<rank>>
    {
      using type = typename ::FEValuesViews::Tensor<rank, dim, spacedim>;
    };
 
    template <int dim, int spacedim, int rank>
    struct ViewType<dim, spacedim, FEValuesExtractors::SymmetricTensor<rank>>
    {
      using type =
        typename ::FEValuesViews::SymmetricTensor<rank, dim, spacedim>;
    };
 
    template <int dim, int spacedim>
    struct Cache
    {
      std::vector<Lazy<::FEValuesViews::Scalar<dim, spacedim>>> scalars;
      std::vector<Lazy<::FEValuesViews::Vector<dim, spacedim>>> vectors;
      std::vector<
        Lazy<::FEValuesViews::SymmetricTensor<2, dim, spacedim>>>
        symmetric_second_order_tensors;
      std::vector<Lazy<::FEValuesViews::Tensor<2, dim, spacedim>>>
        second_order_tensors;
 
      Cache(const FEValuesBase<dim, spacedim> &fe_values);
    };
  } // namespace FEValuesViews
} // namespace internal
 
namespace FEValuesViews
{
  template <int dim, int spacedim, typename Extractor>
  using View = typename ::internal::FEValuesViews::
    ViewType<dim, spacedim, Extractor>::type;
} // namespace FEValuesViews
 
#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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
    AssertIndexRange(shape_function, fe_values->fe->n_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
 
    AssertIndexRange(shape_function, fe_values->fe->n_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
    DEAL_II_ASSERT_UNREACHABLE();
    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
    AssertIndexRange(shape_function, fe_values->fe->n_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
    AssertIndexRange(shape_function, fe_values->fe->n_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 internal
  {
    inline ::SymmetricTensor<2, 1>
    symmetrize_single_row(const unsigned int n, const ::Tensor<1, 1> &t)
    {
      AssertIndexRange(n, 1);
      (void)n;
 
      return {{t[0]}};
    }
 
 
 
    inline ::SymmetricTensor<2, 2>
    symmetrize_single_row(const unsigned int n, const ::Tensor<1, 2> &t)
    {
      switch (n)
        {
          case 0:
            {
              return {{t[0], 0, t[1] / 2}};
            }
          case 1:
            {
              return {{0, t[1], t[0] / 2}};
            }
          default:
            {
              AssertIndexRange(n, 2);
              return {};
            }
        }
    }
 
 
 
    inline ::SymmetricTensor<2, 3>
    symmetrize_single_row(const unsigned int n, const ::Tensor<1, 3> &t)
    {
      switch (n)
        {
          case 0:
            {
              return {{t[0], 0, 0, t[1] / 2, t[2] / 2, 0}};
            }
          case 1:
            {
              return {{0, t[1], 0, t[0] / 2, 0, t[2] / 2}};
            }
          case 2:
            {
              return {{0, 0, t[2], 0, t[0] / 2, t[1] / 2}};
            }
          default:
            {
              AssertIndexRange(n, 3);
              return {};
            }
        }
    }
  } // namespace internal
 
 
 
  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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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 internal::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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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 \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_j}.
        // (which is incidentally also
        // b_j \dealcoloneq \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 \dealcoloneq \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 \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_ii}.
        // b_jj \dealcoloneq \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
      {
        DEAL_II_NOT_IMPLEMENTED();
        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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
      {
        DEAL_II_NOT_IMPLEMENTED();
        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
  {
    AssertIndexRange(shape_function, fe_values->fe->n_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
      {
        DEAL_II_NOT_IMPLEMENTED();
        gradient_type return_value;
        return return_value;
      }
  }
} // namespace FEValuesViews
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif