MacroGridFactory.hpp

#pragma once

#include <dune/common/typeutilities.hh>
#include <dune/common/version.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/yaspgrid.hh>

namespace Dune
{
  template <class GridType>
  struct CubedWrapper : public GridType {};


  template <class GridType>
  class GridFactory<CubedWrapper<GridType>>
      : public GridFactoryInterface<GridType>
  {
    using Self = GridFactory;
    using Super = GridFactoryInterface<GridType>;
    using ctype = typename GridType::ctype;

    enum { dim = GridType::dimension };
    enum { dimworld = GridType::dimensionworld };

  public:
    template <class... Args, disableCopyMove<Self, Args...> = 0>
    GridFactory (Args&&... args)
      : factory_(std::make_shared<GridFactory<GridType>>(std::forward<Args>(args)...))
    {}

    GridFactory (GridFactory<GridType>& factory)
      : factory_(stackobject_to_shared_ptr(factory))
    {}

    void insertVertex (const FieldVector<ctype,dimworld>& pos) override
    {
      factory_->insertVertex(pos);
    }


    void insertElement (const GeometryType& type,
                        const std::vector<unsigned int>& vertices) override
    {
      // triangulation of reference cube
      static const auto reference_cubes = std::make_tuple(
        std::array<std::array<int,2>, 1>{std::array<int,2>{0,1}},
        std::array<std::array<int,3>, 2>{std::array<int,3>{3,0,1}, std::array<int,3>{0,3,2}},
        std::array<std::array<int,4>, 6>{std::array<int,4>{0,7,3,1}, std::array<int,4>{0,7,5,1},
                                         std::array<int,4>{0,7,5,4}, std::array<int,4>{0,7,3,2},
                                         std::array<int,4>{0,7,6,2}, std::array<int,4>{0,7,6,4}} );

      assert(type == GeometryTypes::cube(dim));

      auto const& simplices = std::get<dim-1>(reference_cubes);
      thread_local std::vector<unsigned int> corners(dim+1);
      for (auto const& simplex : simplices) {
        for (std::size_t i = 0; i < simplex.size(); ++i)
          corners[i] = vertices[simplex[i]];

        factory_->insertElement(GeometryTypes::simplex(dim), corners);
      }
    }

    using Super::insertElement;

    // TODO: maybe split boundary segment in simplices
    void insertBoundarySegment (const std::vector<unsigned int>& vertices) override
    {
      factory_->insertBoundarySegment(vertices);
    }

    using Super::insertBoundarySegment;

    GridType* createGrid () override
    {
      return factory_->createGrid();
    }

  private:
    std::shared_ptr<GridFactory<GridType>> factory_;
  };

} // end namespace Dune


namespace AMDiS
{
  template <class GridType>
  class MacroGridFactory
  {
    using ctype = typename GridType::ctype;

    enum { dim = GridType::dimension };
    enum { dimworld = GridType::dimensionworld };

  public:

    static std::unique_ptr<GridType> createSimplexGrid (const Dune::FieldVector<ctype,dimworld>& lowerLeft,
                                                        const Dune::FieldVector<ctype,dimworld>& upperRight,
                                                        const std::array<unsigned int,dim>& numElements)
    {
      Dune::GridFactory<Dune::CubedWrapper<GridType>> factory;

      // fallback implementation using temporary YaspGrid
      using TempGrid = Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<double,dim>>;
      auto grid = Dune::StructuredGridFactory<TempGrid>::createCubeGrid(lowerLeft, upperRight, numElements);
      for (auto const& v : vertices(grid->leafGridView()))
        factory.insertVertex(v.geometry().corner(0));

      auto const& indexSet = grid->leafIndexSet();
      for (auto const& e : elements(grid->leafGridView()))
      {
        thread_local std::vector<unsigned int> vertices;
        vertices.resize(e.subEntities(dim));
        for (unsigned int i = 0; i < e.subEntities(dim); ++i)
          vertices[i] = indexSet.subIndex(e,i,dim);
        factory.insertElement(Dune::GeometryTypes::cube(dim), vertices);
      }
      return std::unique_ptr<GridType>(factory.createGrid());
    }
  };

} // end namespace AMDiS