lagrangegridcreator.hh

#pragma once
 
#include <cassert>
#include <cstdint>
#include <limits>
#include <optional>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/geometry/utility/typefromvertexcount.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/localfunctions/lagrange.hh>
#include <dune/grid/common/gridfactory.hh>
 
#include <dune/vtk/types.hh>
#include <dune/vtk/gridcreatorinterface.hh>
#include <dune/vtk/gridfunctions/lagrangegridfunction.hh>
#include <dune/vtk/utility/lagrangepoints.hh>
 
namespace Dune::Vtk
{
    // \brief Create a grid from data that represents higher (Lagrange) cells.
    template <class GridType>
    class LagrangeGridCreator
        : public GridCreatorInterface<GridType, LagrangeGridCreator<GridType>>
    {
      using Self = LagrangeGridCreator;
      using Super = GridCreatorInterface<GridType, Self>;
 
    public:
      using Grid = GridType;
      using Element = typename Grid::template Codim<0>::Entity;
      using GlobalCoordinate = typename Super::GlobalCoordinate;
      using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
 
    private:
      struct ElementParametrization
      {
        GeometryType type = GeometryTypes::none(Grid::dimension); //< Geometry type of the element
        std::vector<std::int64_t> nodes = {};    //< Indices of the w.r.t. `nodes_` vector
        std::vector<unsigned int> corners = {};  //< Insertion-indices of the element corner nodes
      };
 
      using Parametrization = std::vector<ElementParametrization>;
      using Nodes = std::vector<GlobalCoordinate>;
 
      class LocalParametrization;
      class LocalFunction;
 
    public:
      using LocalGeometry = MultiLinearGeometry<typename Element::Geometry::ctype,Element::dimension,Element::dimension>;
 
    public:
      using Super::Super;
      using Super::factory;
 
      void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
                              std::vector<std::uint64_t> const& /*point_ids*/)
      {
        // store point coordinates in member variable
        nodes_ = points;
      }
 
      template <class F>
      using HasParametrizedElements = decltype(std::declval<F>().insertElement(std::declval<GeometryType>(),
        std::declval<std::vector<unsigned int> const&>(), std::declval<std::function<GlobalCoordinate(LocalCoordinate)>>()));
 
      void insertElementsImpl (std::vector<std::uint8_t> const& types,
                              std::vector<std::int64_t> const& offsets,
                              std::vector<std::int64_t> const& connectivity)
      {
        assert(nodes_.size() > 0);
 
        // mapping of node index to element-vertex index
        std::vector<std::int64_t> elementVertices(nodes_.size(), -1);
        parametrization_.reserve(types.size());
 
        std::int64_t vertexIndex = 0;
        for (std::size_t i = 0; i < types.size(); ++i) {
          auto type = Vtk::to_geometry(types[i]);
          if (type.dim() != Grid::dimension)
            continue;
 
          Vtk::CellType cellType{type};
          auto refElem = referenceElement<double,Grid::dimension>(type);
 
          std::int64_t shift = (i == 0 ? 0 : offsets[i-1]);
          std::int64_t nNodes = offsets[i] - shift;
          int nVertices = refElem.size(Grid::dimension);
 
          // insert vertices into grid and construct element vertices
          std::vector<unsigned int> element(nVertices);
          for (int j = 0; j < nVertices; ++j) {
            auto index = connectivity.at(shift + j);
            auto& vertex = elementVertices.at(index);
            if (vertex < 0) {
              factory().insertVertex(nodes_.at(index));
              vertex = vertexIndex++;
            }
            element[j] = vertex;
          }
 
          // permute element indices
          if (!cellType.noPermutation()) {
            // apply index permutation
            std::vector<unsigned int> cell(element.size());
            for (int j = 0; j < int(element.size()); ++j)
              cell[j] = element[cellType.permutation(j)];
            std::swap(element, cell);
          }
 
          // fill vector of element parametrizations
          parametrization_.push_back(ElementParametrization{type});
          auto& param = parametrization_.back();
 
          param.nodes.resize(nNodes);
          for (std::int64_t j = 0; j < nNodes; ++j)
            param.nodes[j] = connectivity.at(shift + j);
          param.corners = element;
 
          // try to create element with parametrization
          if constexpr (Std::is_detected_v<HasParametrizedElements, GridFactory<Grid>>) {
            try {
              factory().insertElement(type, element,
                localParametrization(parametrization_.size()-1));
            } catch (Dune::GridError const& /* notImplemented */) {
              factory().insertElement(type, element);
            }
          } else {
            factory().insertElement(type, element);
          }
        }
      }
 
 
      LocalParametrization localParametrization (unsigned int insertionIndex) const
      {
        VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
        auto const& localParam = parametrization_.at(insertionIndex);
        return LocalParametrization{nodes_, localParam, order(localParam)};
      }
 
 
      LocalParametrization localParametrization (Element const& element) const
      {
        VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
 
        unsigned int insertionIndex = factory().insertionIndex(element);
        auto const& localParam = parametrization_.at(insertionIndex);
        VTK_ASSERT(element.type() == localParam.type);
 
        return {nodes_, localParam, order(localParam), localGeometry(element, localParam)};
      }
 
 
      LocalGeometry localGeometry (Element const& element) const
      {
        VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
 
        unsigned int insertionIndex = factory().insertionIndex(element);
        auto const& localParam = parametrization_.at(insertionIndex);
        VTK_ASSERT(element.type() == localParam.type);
 
        return localGeometry(element, localParam);
      }
 
    private:
      // implementation details of localGeometry()
      LocalGeometry localGeometry (Element const& element, ElementParametrization const& localParam) const
      {
        auto refElem = referenceElement<typename Element::Geometry::ctype,Element::dimension>(localParam.type);
 
        // collect indices of vertices
        std::vector<unsigned int> indices(refElem.size(Grid::dimension));
        for (int i = 0; i < refElem.size(Grid::dimension); ++i)
          indices[i] = factory().insertionIndex(element.template subEntity<Grid::dimension>(i));
 
        // calculate permutation vector
        std::vector<unsigned int> permutation(indices.size());
        for (std::size_t i = 0; i < indices.size(); ++i) {
          auto it = std::find(localParam.corners.begin(), localParam.corners.end(), indices[i]);
          VTK_ASSERT(it != localParam.corners.end());
          permutation[i] = std::distance(localParam.corners.begin(), it);
        }
 
        std::vector<LocalCoordinate> corners(permutation.size());
        for (std::size_t i = 0; i < permutation.size(); ++i)
          corners[i] = refElem.position(permutation[i], Element::dimension);
 
        return {localParam.type, corners};
      }
 
    public:
 
      int order (GeometryType type, std::size_t nNodes) const
      {
        for (int o = 1; o <= int(nNodes); ++o)
          if (numLagrangePoints(type, o) == std::size_t(nNodes))
            return o;
 
        return 1;
      }
 
      int order (ElementParametrization const& localParam) const
      {
        return order(localParam.type, localParam.nodes.size());
      }
 
      int order () const
      {
        VTK_ASSERT(!parametrization_.empty());
        auto const& localParam = parametrization_.front();
        return order(localParam);
      }
 
    public:
 
      friend LocalFunction localFunction (LagrangeGridCreator& gridCreator)
      {
        return LocalFunction{gridCreator};
      }
 
      friend LocalFunction localFunction (LagrangeGridCreator const& gridCreator)
      {
        return LocalFunction{gridCreator};
      }
 
      friend LocalFunction localFunction (LagrangeGridCreator&& gridCreator)
      {
        DUNE_THROW(Dune::Exception, "Cannot pass temporary LagrangeGridCreator to localFunction(). Pass an lvalue-reference instead.");
        return LocalFunction{gridCreator};
      }
 
      std::string name () const
      {
        return "LagrangeParametrization";
      }
 
      int numComponents () const
      {
        return GlobalCoordinate::size();
      }
 
      Vtk::DataTypes dataType () const
      {
        return dataTypeOf<typename Element::Geometry::ctype>();
      }
 
      struct EntitySet
      {
        using Grid = GridType;
        using GlobalCoordinate = typename Self::GlobalCoordinate;
      };
 
      EntitySet entitySet () const
      {
        assert(false && "Should not be used!");
        return EntitySet{};
      }
 
      GlobalCoordinate operator() (GlobalCoordinate const&) const
      {
        assert(false && "Should not be used!");
        return GlobalCoordinate{};
      }
 
    private:
      Nodes nodes_;
 
      Parametrization parametrization_;
    };
 
    // deduction guides
    template <class Grid>
    LagrangeGridCreator(GridFactory<Grid>&)
      -> LagrangeGridCreator<Grid>;
 
    template <class GridType, class Range, class Context>
    struct AssociatedGridFunction<LagrangeGridCreator<GridType>, Range, Context>
    {
      using type = LagrangeGridFunction<GridType, Range, Context>;
    };
 
    template <class Grid>
    class LagrangeGridCreator<Grid>::LocalParametrization
    {
      using ctype = typename Grid::ctype;
 
      using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
      using LocalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::LocalCoordinate;
      using LocalGeometry = MultiLinearGeometry<ctype,Grid::dimension,Grid::dimension>;
 
      using LocalFE = LagrangeLocalFiniteElement<Vtk::LagrangePointSet, Grid::dimension, ctype, ctype>;
      using LocalBasis = typename LocalFE::Traits::LocalBasisType;
      using LocalBasisTraits = typename LocalBasis::Traits;
 
    public:
      template <class Nodes, class LocalParam>
      LocalParametrization (Nodes const& nodes, LocalParam const& param, int order)
        : localFE_(param.type, order)
        , localNodes_(param.nodes.size())
      {
        for (std::size_t i = 0; i < localNodes_.size(); ++i)
          localNodes_[i] = nodes[param.nodes[i]];
      }
 
      template <class Nodes, class LocalParam, class LG>
      LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, LG&& localGeometry)
        : LocalParametrization(nodes, param, order)
      {
        localGeometry_.emplace(std::forward<LG>(localGeometry));
      }
 
      template <class LocalCoordinate>
      GlobalCoordinate operator() (LocalCoordinate const& local) const
      {
        // map coordinates if element corners are permuted
        LocalCoordinate x = localGeometry_ ? localGeometry_->global(local) : local;
 
        LocalBasis const& localBasis = localFE_.localBasis();
        localBasis.evaluateFunction(x, shapeValues_);
        assert(shapeValues_.size() == localNodes_.size());
 
        using field_type = typename LocalBasisTraits::RangeType::field_type;
 
        GlobalCoordinate out(0);
        for (std::size_t i = 0; i < shapeValues_.size(); ++i)
          out.axpy(field_type(shapeValues_[i]), localNodes_[i]);
 
        return out;
      }
 
    private:
      LocalFE localFE_;
      std::vector<GlobalCoordinate> localNodes_;
      std::optional<LocalGeometry> localGeometry_;
 
      mutable std::vector<typename LocalBasisTraits::RangeType> shapeValues_;
    };
 
 
    template <class Grid>
    class LagrangeGridCreator<Grid>::LocalFunction
    {
      using ctype = typename Grid::ctype;
      using LocalContext = typename Grid::template Codim<0>::Entity;
      using GlobalCoordinate = typename LocalContext::Geometry::GlobalCoordinate;
      using LocalCoordinate = typename LocalContext::Geometry::LocalCoordinate;
      using LocalParametrization = typename LagrangeGridCreator::LocalParametrization;
 
    public:
      explicit LocalFunction (LagrangeGridCreator const& gridCreator)
        : gridCreator_(&gridCreator)
      {}
 
      explicit LocalFunction (LagrangeGridCreator&& gridCreator) = delete;
 
      void bind (LocalContext const& element)
      {
        localContext_ = element;
        localParametrization_.emplace(gridCreator_->localParametrization(element));
      }
 
      void unbind () { /* do nothing */ }
 
      GlobalCoordinate operator() (LocalCoordinate const& local) const
      {
        assert(!!localParametrization_);
        return (*localParametrization_)(local);
      }
 
      LocalContext const& localContext () const
      {
        return localContext_;
      }
 
    private:
      LagrangeGridCreator const* gridCreator_;
 
      LocalContext localContext_;
      std::optional<LocalParametrization> localParametrization_;
    };
 
} // end namespace Dune::Vtk