fe_field_function.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2007 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_fe_function_h
#define dealii_fe_function_h

#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/thread_local_storage.h>

#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/mapping_q1.h>

#include <deal.II/grid/grid_tools_cache.h>

#include <deal.II/lac/vector.h>

#include <boost/optional.hpp>


DEAL_II_NAMESPACE_OPEN

namespace VectorTools
{
  class ExcPointNotAvailableHere;
}

namespace Functions
{
  template <int dim,
            typename DoFHandlerType = DoFHandler<dim>,
            typename VectorType     = Vector<double>>
  class FEFieldFunction : public Function<dim, typename VectorType::value_type>
  {
  public:
    FEFieldFunction(
      const DoFHandlerType &dh,
      const VectorType &    data_vector,
      const Mapping<dim> &  mapping = StaticMappingQ1<dim>::mapping);

    void
    set_active_cell(
      const typename DoFHandlerType::active_cell_iterator &newcell);

    virtual void
    vector_value(
      const Point<dim> &                       p,
      Vector<typename VectorType::value_type> &values) const override;

    virtual typename VectorType::value_type
    value(const Point<dim> &p, const unsigned int component = 0) const override;

    virtual void
    value_list(const std::vector<Point<dim>> &               points,
               std::vector<typename VectorType::value_type> &values,
               const unsigned int component = 0) const override;


    virtual void
    vector_value_list(const std::vector<Point<dim>> &points,
                      std::vector<Vector<typename VectorType::value_type>>
                        &values) const override;

    virtual void
    vector_gradient(const Point<dim> &p,
                    std::vector<Tensor<1, dim, typename VectorType::value_type>>
                      &gradients) const override;

    virtual Tensor<1, dim, typename VectorType::value_type>
    gradient(const Point<dim> & p,
             const unsigned int component = 0) const override;

    virtual void
    vector_gradient_list(
      const std::vector<Point<dim>> &p,
      std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
        &gradients) const override;

    virtual void
    gradient_list(
      const std::vector<Point<dim>> &                               p,
      std::vector<Tensor<1, dim, typename VectorType::value_type>> &gradients,
      const unsigned int component = 0) const override;


    virtual typename VectorType::value_type
    laplacian(const Point<dim> & p,
              const unsigned int component = 0) const override;

    virtual void
    vector_laplacian(
      const Point<dim> &                       p,
      Vector<typename VectorType::value_type> &values) const override;

    virtual void
    laplacian_list(const std::vector<Point<dim>> &               points,
                   std::vector<typename VectorType::value_type> &values,
                   const unsigned int component = 0) const override;

    virtual void
    vector_laplacian_list(const std::vector<Point<dim>> &points,
                          std::vector<Vector<typename VectorType::value_type>>
                            &values) const override;

    unsigned int
    compute_point_locations(
      const std::vector<Point<dim>> &                             points,
      std::vector<typename DoFHandlerType::active_cell_iterator> &cells,
      std::vector<std::vector<Point<dim>>> &                      qpoints,
      std::vector<std::vector<unsigned int>> &                    maps) const;

  private:
    using cell_hint_t = Threads::ThreadLocalStorage<
      typename DoFHandlerType::active_cell_iterator>;

    SmartPointer<const DoFHandlerType,
                 FEFieldFunction<dim, DoFHandlerType, VectorType>>
      dh;

    const VectorType &data_vector;

    const Mapping<dim> &mapping;

    GridTools::Cache<dim, DoFHandlerType::space_dimension> cache;

    mutable cell_hint_t cell_hint;

    boost::optional<Point<dim>>
    get_reference_coordinates(
      const typename DoFHandlerType::active_cell_iterator &cell,
      const Point<dim> &                                   point) const;
  };
} // namespace Functions

DEAL_II_NAMESPACE_CLOSE

#endif