gridfunction.hh

#pragma once
 
#include <memory>
#include <numeric>
#include <string>
#include <type_traits>
#include <vector>
 
#include <dune/common/concept.hh>
#include <dune/common/referencehelper.hh>
#include <dune/common/typetraits.hh>
#include <dune/vtk/types.hh>
#include <dune/vtk/gridfunctions/componentmapper.hh>
#include <dune/vtk/gridfunctions/localfunction.hh>
#include <dune/vtk/gridfunctions/rangeproxy.hh>
#include <dune/vtk/gridfunctions/rangetype.hh>
#include <dune/vtk/gridfunctions/vtkfunctionwrapper.hh>
 
namespace Dune::Vtk {
namespace Impl {
 
template <class GridViewType, class ComponentMapperType>
struct GridFunctionDefinition
{
  using GridView = GridViewType;
  using ComponentMapper = ComponentMapperType;
  using Element = typename GridView::template Codim<0>::Entity;
  using Domain = typename Element::Geometry::GlobalCoordinate;
  using LocalFunction = Vtk::LocalFunction<Element,ComponentMapper>;
 
  // Definition of the interface a LocalFunction must fulfill
  struct Interface
  {
    virtual ~Interface() = default;
    virtual Interface* clone () const = 0;
    virtual double evaluate (int, Domain const&) const = 0;
    virtual LocalFunction asLocalFunction(std::vector<int> const&) const = 0;
 
  };
 
  // Templatized implementation of the interface
  template <class Impl>
  struct Model : public Interface
  {
    template <class I, disableCopyMove<Model,I> = 0>
    Model(I&& impl) : impl_(std::forward<I>(impl)) {}
    Model* clone() const final { return new Model(*this); }
    auto& impl() { return resolveRef(impl_); }
    const auto& impl() const { return resolveRef(impl_); }
 
    double evaluate(int i, Domain const& x) const final {
      return ComponentMapper::map(i, x, impl()); }
    LocalFunction asLocalFunction(std::vector<int> const& comp) const final {
      return LocalFunction{localFunction(impl()), comp};
    }
 
  private:
    Impl impl_;
  };
 
  // The Concept definition of a LocalFunction
  struct Concept
  {
    using LFConcept = typename LocalFunctionDefinition<Element,ComponentMapper>::Concept;
 
    template <class F>
    auto require(F&& f) -> decltype
    (
      Dune::Concept::requireConcept<LFConcept>(localFunction(f))
    );
  };
};
 
} // end namespace Impl
 
 
 
 
template <class GridView, class ComponentMapper = DefaultComponentMapper>
class GridFunction
{
  using Definition = Impl::GridFunctionDefinition<GridView,ComponentMapper>;
  using Interface = typename Definition::Interface;
  using Domain = typename Definition::Domain;
  using LocalFunction = typename Definition::LocalFunction;
 
  template <class F>
  using Model = typename Definition::template Model<std::remove_reference_t<F>>;
 
public:
  template <class F,
    class R = RangeType<ResolveRef_t<F>,Domain>>
  GridFunction (F&& f, std::string name, int numComponents = dimRange<R>())
    : GridFunction(std::forward<F>(f), Vtk::FieldInfo{std::move(name), numComponents, RangeTypes::UNSPECIFIED, DataTypes::FLOAT64})
  {}
 
  template <class F, disableCopyMove<GridFunction, F> = 0,
    decltype(std::declval<F>().name(), bool{}) = true,
    decltype(std::declval<F>().numComponents(), bool{}) = true,
    decltype(std::declval<F>().dataType(), bool{}) = true>
  GridFunction (F&& f)
    : GridFunction(std::forward<F>(f), Vtk::FieldInfo{f.name(), f.numComponents(), Vtk::RangeTypes::UNSPECIFIED, f.dataType()})
  {}
 
  template <class F>
  GridFunction (F&& f, const Vtk::FieldInfo& info)
    : interface_(std::make_unique<Model<F>>(std::forward<F>(f)))
    , name_(info.name())
  {
    static_assert(models<typename Definition::Concept,ResolveRef_t<F>>(),
      "Implementation does not model the Function concept.");
 
    setComponents(info.size());
    setRangeType(info.rangeType());
    setDataType(info.dataType());
  }
 
  GridFunction () = default;
  GridFunction (GridFunction &&) = default;
  GridFunction& operator= (GridFunction&&) = default;
 
  GridFunction (const GridFunction& other)
    : interface_(other.interface_ ? other.interface_->clone() : nullptr)
    , name_(other.name_)
    , components_(other.components_)
    , dataType_(other.dataType_)
    , rangeType_(other.rangeType_)
  {}
 
  GridFunction& operator= (const GridFunction& other)
  {
    return *this = GridFunction(other);
  }
 
  Impl::RangeProxy<GridFunction,Domain> operator() (Domain const& x) const
  {
    return {this, x};
  }
 
  double evaluate (int i, Domain const& x) const
  {
    return i < int(components_.size()) ? interface_->evaluate(components_[i], x) : 0.0;
  }
 
  friend LocalFunction localFunction (GridFunction const& self)
  {
    return self.interface_->asLocalFunction(self.components_);
  }
 
  std::string const& name () const
  {
    return name_;
  }
 
  void setName (std::string name)
  {
    name_ = std::move(name);
  }
 
  int numComponents () const
  {
    return rangeType_ == Vtk::RangeTypes::SCALAR ? 1 :
           rangeType_ == Vtk::RangeTypes::VECTOR ? 3 :
           rangeType_ == Vtk::RangeTypes::TENSOR ? 9 : int(components_.size());
  }
 
  const std::vector<int>& components () const
  {
    return components_;
  }
 
  void setComponents (std::vector<int> components)
  {
    components_ = components;
  }
 
  void setComponents (int ncomps)
  {
    setComponents(allComponents(ncomps));
  }
 
  Vtk::DataTypes dataType () const
  {
    return dataType_;
  }
 
  void setDataType (Vtk::DataTypes type)
  {
    dataType_ = type;
  }
 
  Vtk::RangeTypes rangeType () const
  {
    return rangeType_;
  }
 
  void setRangeType (Vtk::RangeTypes type, std::size_t ncomp = 1)
  {
    rangeType_ = type;
    if (type == Vtk::RangeTypes::AUTO)
      rangeType_ = rangeTypeOf(ncomp);
  }
 
  void setFieldInfo (const Vtk::FieldInfo& info)
  {
    setName(info.name());
    setComponents(info.size());
    setRangeType(info.rangeType());
    setDataType(info.dataType());
  }
 
private:
  std::vector<int> allComponents (int n) const
  {
    std::vector<int> components(n);
    std::iota(components.begin(), components.end(), 0);
    return components;
  }
 
private:
  std::unique_ptr<Interface> interface_;
  std::string name_{};
  std::vector<int> components_{};
  Vtk::DataTypes dataType_ = Vtk::DataTypes::FLOAT64;
  Vtk::RangeTypes rangeType_ = Vtk::RangeTypes::UNSPECIFIED;
};
 
} // end namespace Dune::Vtk