gridfunction.hh

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <array>
#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(), rangeTypeImpl(f, Dune::PriorityTag<1>{}), 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())
    , shape_(staticShape<ResolveRef_t<F>>())
  {
    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_)
    , shape_(other.shape_)
    , outputMap_(other.outputMap_)
    , 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(outputMap_.size()) && outputMap_[i] >= 0
      ? interface_->evaluate(outputMap_[i], x) : 0.0;
  }
 
  friend LocalFunction localFunction (GridFunction const& self)
  {
    return self.interface_->asLocalFunction(self.outputMap_);
  }
 
  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;
    updateOutputMap();
  }
 
  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);
    updateOutputMap();
  }
 
  void setFieldInfo (const Vtk::FieldInfo& info)
  {
    setName(info.name());
    setComponents(info.size());
    setRangeType(info.rangeType());
    setDataType(info.dataType());
  }
 
private:
  template <class F>
  static constexpr std::array<int,2> staticShape ()
  {
    if constexpr (std::is_invocable_v<F const&, Domain const&>)
      return rangeShape<RangeType<F const&,Domain>>();
    else
      return {0,0};
  }
 
  // Shape {rows,cols} of the (selected) range interpreted as a tensor
  std::array<int,2> tensorShape () const
  {
    const int n = components_.size();
    if (shape_[0] > 0 && shape_[0]*shape_[1] <= n)
      return shape_;
    for (int k : {1,2,3})
      if (k*k == n) return {k,k};
    return {(n+2)/3, 3}; // unknown: fill the 3x3 tensor row-wise
  }
 
  // Map each output slot [0,numComponents()) to a range component, or -1 for zero padding
  void updateOutputMap ()
  {
    outputMap_.assign(numComponents(), -1);
    if (rangeType_ == Vtk::RangeTypes::TENSOR) {
      auto [rows,cols] = tensorShape();
      for (int r = 0; r < std::min(rows,3); ++r)
        for (int c = 0; c < std::min(cols,3); ++c)
          if (int k = r*cols + c; k < int(components_.size()))
            outputMap_[3*r + c] = components_[k];
    } else {
      for (std::size_t i = 0; i < std::min(outputMap_.size(), components_.size()); ++i)
        outputMap_[i] = components_[i];
    }
  }
 
  std::vector<int> allComponents (int n) const
  {
    std::vector<int> components(n);
    std::iota(components.begin(), components.end(), 0);
    return components;
  }
 
  // Forward the range-type of `f` if available
  template <class F>
  static auto rangeTypeImpl (F const& f, Dune::PriorityTag<1>)
    -> decltype(Vtk::RangeTypes(f.rangeType()))
  {
    return f.rangeType();
  }
 
  template <class F>
  static Vtk::RangeTypes rangeTypeImpl (F const& /*f*/, Dune::PriorityTag<0>)
  {
    return Vtk::RangeTypes::UNSPECIFIED;
  }
 
private:
  std::unique_ptr<Interface> interface_;
  std::string name_{};
  std::vector<int> components_{};
  std::array<int,2> shape_{0,0};
  std::vector<int> outputMap_{};
  Vtk::DataTypes dataType_ = Vtk::DataTypes::FLOAT64;
  Vtk::RangeTypes rangeType_ = Vtk::RangeTypes::UNSPECIFIED;
};
 
} // end namespace Dune::Vtk