lagrangegridfunction.hh

Go to the documentation of this file.

#pragma once
 
#include <optional>
#include <vector>
 
#include <dune/common/dynvector.hh>
#include <dune/common/version.hh>
#include <dune/localfunctions/lagrange.hh>
#include <dune/vtk/gridfunctions/common.hh>
#include <dune/vtk/gridfunctions/continuousgridfunction.hh>
#include <dune/vtk/utility/errors.hh>
#include <dune/vtk/utility/lagrangepoints.hh>
 
namespace Dune::Vtk
{
  template <class GridType>
  class LagrangeGridCreator;
 
  template <class GridType, class RangeType, class Context>
  class LagrangeGridFunction
  {
    using Grid = GridType;
    using Field = typename Dune::FieldTraits<RangeType>::field_type;
 
    using Factory = GridFactory<Grid>;
    using GridCreator = LagrangeGridCreator<Grid>;
 
    // The range must fulfill the DenseVector interface. Scalar types
    // are simply wrapped into a FieldVector<T,1> on extraction. This
    // can then be converted into the type `T` on evaluation.
    using RangeWrapperType = std::conditional_t<std::is_arithmetic_v<RangeType>,
      Dune::FieldVector<RangeType,1>, RangeType>;
 
  public:
    struct EntitySet
    {
      using Grid = GridType;
      using Element = typename GridType::template Codim<0>::Entity;
      using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
      using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
    };
 
    using Domain = typename EntitySet::GlobalCoordinate;
    using Range = RangeType;
    using Signature = Range(Domain);
 
  private:
    template <class LC>
    class PointDataLocalFunction
    {
      using LFE = LagrangeLocalFiniteElement<LagrangePointSet, LC::dimension, Field, Field>;
      using LB = typename LFE::Traits::LocalBasisType;
 
    public:
      using LocalContext = LC;
      using Domain = typename LC::Geometry::LocalCoordinate;
      using Range = RangeType;
      using Signature = Range(Domain);
 
    public:
      PointDataLocalFunction (GridCreator const* creator, std::vector<Field> const* values, unsigned int comp,
                              std::vector<std::uint8_t> const* types,
                              std::vector<std::int64_t> const* offsets,
                              std::vector<std::int64_t> const* connectivity)
        : creator_(creator)
        , values_(values)
        , comp_(comp)
        , types_(types)
        , offsets_(offsets)
        , connectivity_(connectivity)
      {}
 
      PointDataLocalFunction () = default;
 
 
      void bind (LocalContext const& element)
      {
        unsigned int insertionIndex = creator_->factory().insertionIndex(element);
 
        std::int64_t shift = (insertionIndex == 0 ? 0 : (*offsets_)[insertionIndex-1]);
        std::int64_t numNodes = (*offsets_)[insertionIndex] - shift;
        [[maybe_unused]] std::int64_t maxNumNodes = numLagrangePoints(element.type(), 20);
        VTK_ASSERT(numNodes > 0 && numNodes < maxNumNodes);
 
        int order = creator_->order(element.type(), numNodes);
        VTK_ASSERT(order > 0 && order < 20);
 
        // construct a local finite-element with the corresponding order and Lagrange points
        // as stored in the file
        lfe_.emplace(LFE{element.type(), (unsigned int)(order)});
        lgeo_.emplace(creator_->localGeometry(element));
 
        // collect values on lagrange nodes
        localValues_.resize(numNodes);
        for (std::int64_t i = shift, i0 = 0; i < (*offsets_)[insertionIndex]; ++i, ++i0) {
          std::int64_t idx = (*connectivity_)[i];
          auto& v = localValues_[i0];
          Impl::resizeRange(v,comp_);
          for (unsigned int j = 0; j < comp_; ++j)
            v[j] = (*values_)[comp_*idx + j];
        }
      }
 
      void unbind ()
      {
        lfe_.reset();
        lgeo_.reset();
      }
 
      // NOTE: do we need to transform the local coordinates?
      Range operator() (Domain const& local) const
      {
        assert(!!lfe_);
        auto const& lb = lfe_->localBasis();
        lb.evaluateFunction(lgeo_->global(local), shapeValues_);
        assert(shapeValues_.size() == localValues_.size());
 
        auto y = Impl::constructRange<RangeWrapperType>(comp_, 0);
        for (std::size_t i = 0; i < shapeValues_.size(); ++i)
          y.axpy(shapeValues_[i], localValues_[i]);
 
        return Range(y);
      }
 
      friend DummyLocalFunctionDerivative<LocalContext,Range,Domain> derivative (PointDataLocalFunction const&) { return {}; }
 
    private:
      GridCreator const* creator_ = nullptr;
      std::vector<Field> const* values_ = nullptr;
      unsigned int comp_;
      std::vector<std::uint8_t> const* types_ = nullptr;
      std::vector<std::int64_t> const* offsets_ = nullptr;
      std::vector<std::int64_t> const* connectivity_ = nullptr;
 
      // Local Finite-Element
      std::optional<LFE> lfe_ = std::nullopt;
      std::optional<typename GridCreator::LocalGeometry> lgeo_ = std::nullopt;
 
      // cache of local values
      std::vector<RangeWrapperType> localValues_;
      mutable std::vector<typename LB::Traits::RangeType> shapeValues_;
    };
 
    template <class LC>
    class CellDataLocalFunction
    {
    public:
      using LocalContext = LC;
      using Domain = typename LC::Geometry::LocalCoordinate;
      using Range = RangeType;
      using Signature = Range(Domain);
 
    public:
      CellDataLocalFunction (GridCreator const* creator, std::vector<Field> const* values, unsigned int comp,
                              std::vector<std::uint8_t> const* /*types*/,
                              std::vector<std::int64_t> const* /*offsets*/,
                              std::vector<std::int64_t> const* /*connectivity*/)
        : creator_(creator)
        , values_(values)
        , comp_(comp)
      {}
 
      CellDataLocalFunction () = default;
 
      void bind (LocalContext const& element)
      {
        unsigned int idx = creator_->factory().insertionIndex(element);
 
        // collect values on cells
        auto& v = localValue_;
        Impl::resizeRange(v,comp_);
 
        for (unsigned int j = 0; j < comp_; ++j)
          v[j] = (*values_)[comp_*idx + j];
      }
 
      void unbind ()
      {}
 
      Range operator() (Domain const& local) const
      {
        return Range(localValue_);
      }
 
    private:
      GridCreator const* creator_ = nullptr;
      std::vector<Field> const* values_ = nullptr;
      unsigned int comp_;
 
      // cache of local values
      RangeWrapperType localValue_;
    };
 
    template <class LC>
    using LocalFunction = std::conditional_t< std::is_same_v<Context,PointContext>,
      PointDataLocalFunction<LC>,
      CellDataLocalFunction<LC>>;
 
  public:
    LagrangeGridFunction (GridCreator const& creator, std::vector<Field> const& values,
                          std::string name ,unsigned int ncomps, Vtk::DataTypes dataType,
                          std::vector<std::uint8_t> const& types,
                          std::vector<std::int64_t> const& offsets,
                          std::vector<std::int64_t> const& connectivity)
      : creator_(&creator)
      , values_(&values)
      , name_(std::move(name))
      , ncomps_(ncomps)
      , dataType_(dataType)
      , types_(&types)
      , offsets_(&offsets)
      , connectivity_(&connectivity)
    {}
 
    LagrangeGridFunction () = default;
 
    Range operator() (Domain const& global) const
    {
      DUNE_THROW(Dune::NotImplemented, "Evaluation in global coordinates not implemented.");
      return Impl::constructRange<Range>(ncomps_, 0);
    }
 
    EntitySet const& entitySet () const
    {
      return entitySet_;
    }
 
    std::string const& name () const
    {
      return name_;
    }
 
    int numComponents () const
    {
      return ncomps_;
    }
 
    Vtk::DataTypes dataType () const
    {
      return dataType_;
    }
 
    friend LocalFunction<typename EntitySet::Element> localFunction (LagrangeGridFunction const& gf)
    {
      return {gf.creator_, gf.values_, gf.ncomps_, gf.types_, gf.offsets_, gf.connectivity_};
    }
 
  private:
    GridCreator const* creator_ = nullptr;
    std::vector<Field> const* values_ = nullptr;
    std::string name_ = "GridFunction";
    unsigned int ncomps_ = 0;
    Vtk::DataTypes dataType_ = Vtk::DataTypes::UNKNOWN;
    std::vector<std::uint8_t> const* types_ = nullptr;
    std::vector<std::int64_t> const* offsets_ = nullptr;
    std::vector<std::int64_t> const* connectivity_ = nullptr;
 
    EntitySet entitySet_;
  };
 
} // end namespace Dune::Vtk