MeshCreator.hpp

#pragma once

#include <algorithm>
#include <array>
#include <memory>
#include <string>

#include <dune/common/filledarray.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typeutilities.hh>

#if HAVE_ALBERTA
#include <dune/grid/albertagrid/albertareader.hh>
#endif
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/utility/structuredgridfactory.hh>

#include <amdis/AdaptiveGrid.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/Filesystem.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/utility/MacroGridFactory.hpp>

namespace Dune
{
  // forward declarations
  template <int dim, class Coordinates> class YaspGrid;
  template <class ct, int dim> class EquidistantCoordinates;
  template <class ct, int dim> class EquidistantOffsetCoordinates;
  template <class ct, int dim> class TensorProductCoordinates;
}


namespace AMDiS
{
  /// A creator class for dune grids.
  template <class Grid>
  struct MeshCreator
  {
    enum { dimension = Grid::dimension };
    enum { dimworld = Grid::dimensionworld };

    using ctype = typename Grid::ctype;
    using HostGrid = typename AdaptiveGrid_t<Grid>::HostGrid;

    /// Construct a new MeshCreator
    /**
     * \param name  The name of the mesh used in the initfile
     **/
    MeshCreator(std::string const& name)
      : name_(name)
    {}

    /// Create a new grid by inspecting the intifile parameter group `[meshName]`
    /**
     * Reads first the parameter `[meshName]->macro file name` and if set
     * - reads the grid from file
     *
     * Otherwise reads the parameter `[meshName]->structured` and if set, creates:
     * - cube: a structured cube grid
     * - simplex: a structured simplex grid
     * using a StructuredGridFactory
     *
     * Otherwise tries to create a grid depending on the grid type.
     **/
    std::shared_ptr<Grid> create() const
    {
      auto filename = Parameters::get<std::string>(name_ + "->macro file name");
      auto structured = Parameters::get<std::string>(name_ + "->structured");
      std::unique_ptr<HostGrid> gridPtr;
      if (filename) {
        // read a macro file
        gridPtr = create_unstructured_grid(filename.value());
      } else if (structured) {
        // create structured grids
        if (structured.value() == "cube") {
          gridPtr = create_cube_grid();
        } else if (structured && structured.value() == "simplex") {
          gridPtr = create_simplex_grid();
        } else {
          error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
        }
      } else {
        // decide by inspecting the grid type how to create the grid
        gridPtr = create_by_gridtype<HostGrid>(Dune::PriorityTag<42>{});
      }

      // Perform initial refinement and load balance if requested in the initfile
      auto globalRefinements = Parameters::get<int>(name_ + "->global refinements");
      if (globalRefinements)
        gridPtr->globalRefine(globalRefinements.value());
      auto loadBalance = Parameters::get<bool>(name_ + "->load balance");
      if (loadBalance && loadBalance.value())
        gridPtr->loadBalance();

      return construct(std::move(gridPtr));
    }

    /// Create a structured cube grid
    std::unique_ptr<HostGrid> create_cube_grid() const
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
        return Dune::StructuredGridFactory<HostGrid>::createCubeGrid(lower, upper, numCells);
      });
    }

    /// Create a structured simplex grid
    std::unique_ptr<HostGrid> create_simplex_grid() const
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
        return Dune::StructuredGridFactory<HostGrid>::createSimplexGrid(lower, upper, numCells);
      });
    }

    /// Return the filled vector of boundary ids. NOTE: not all creators support reading this.
    std::vector<int> const& boundaryIds() const
    {
      return boundaryIds_;
    }

    /// Return the filled vector of element ids. NOTE: not all creators support reading this.
    std::vector<int> const& elementIds() const
    {
      return elementIds_;
    }

  private:
    static std::shared_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid)
    {
      return std::move(hostGrid);
    }

    template <class HG>
    static std::shared_ptr<Grid> construct(std::unique_ptr<HG> hostGrid)
    {
      return std::make_shared<Grid>(std::move(hostGrid));
    }

    // use the structured grid factory to create the grid
    template <class Size = unsigned int, class Factory>
    std::unique_ptr<HostGrid> create_structured_grid(Factory factory) const
    {
      // Lower left corner of the domain
      Dune::FieldVector<ctype,int(dimworld)> lower(0);
      // Upper right corner of the domain
      Dune::FieldVector<ctype,int(dimworld)> upper(1);
      // number of blocks in each coordinate direction
      auto numCells = Dune::filledArray<std::size_t(dimension),Size>(1);

      Parameters::get(name_ + "->min corner", lower);
      Parameters::get(name_ + "->max corner", upper);
      Parameters::get(name_ + "->num cells", numCells);
      return factory(lower, upper, numCells);
    }

    // read a filename from `[meshName]->macro file name` and determine from the extension the fileformat
    std::unique_ptr<HostGrid> create_unstructured_grid(std::string const& filename) const
    {
      filesystem::path fn(filename);
      auto ext = fn.extension();

      if (ext == ".msh") {
        return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{});
      }
      else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
        return read_alberta_file<HostGrid>(filename, Dune::PriorityTag<42>{});
      }
      else {
        error_exit("Cannot read grid file. Unknown file extension.");
        return {};
      }
    }

    template <class GridType>
    using SupportsGmshReader = decltype(std::declval<Dune::GridFactory<GridType>>().insertBoundarySegment(
      std::declval<std::vector<unsigned int>>(),
      std::declval<std::shared_ptr<Dune::BoundarySegment<GridType::dimension, GridType::dimensionworld> >>()) );

    // use GmshReader if GridFactory supports insertBoundarySegments
    template <class GridType = HostGrid,
      REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)>
    std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<1>) const
    {
      Dune::GmshReader<GridType> reader;
      return std::unique_ptr<GridType>{reader.read(filename)}; // , boundaryIds_, elementIds_)};
    }

    // fallback if GmshReader cannot be used
    template <class GridType = HostGrid>
    std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<0>) const
    {
      error_exit("Gmsh reader not supported for this grid type!");
      return {};
    }

    // read from Alberta file

#if HAVE_ALBERTA
    template <class GridType>
    using IsAlbertaGrid = decltype(std::declval<GridType>().alberta2dune(0,0));

    // construct the albertagrid directly from a filename
    template <class GridType = HostGrid,
      REQUIRES(GridType::dimensionworld == DIM_OF_WORLD),
      REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<3>) const
    {
      return std::make_unique<GridType>(filename);
    }

    // use a gridfactory and the generic AlbertaReader
    template <class GridType = HostGrid,
      REQUIRES(GridType::dimensionworld == DIM_OF_WORLD)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<2>) const
    {
      Dune::GridFactory<GridType> factory;
      if (Environment::mpiRank() == 0) {
        Dune::AlbertaReader<GridType> reader;
        reader.readGrid(filename, factory);
      }
      return std::unique_ptr<GridType>{factory.createGrid()};
    }

    // error if WORLDDIM not the same as Grid::dimensionworld
    template <class GridType = HostGrid,
      REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
    {
      error_exit("AlbertaGrid (and AlbertaReader) require WORLDDIM == Grid::dimensionworld. Change the cmake parameters!");
      return {};
    }
#else
    // fallback if no Alberta is available
    template <class GridType>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<0>) const
    {
      error_exit("AlbertaGrid (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
      return {};
    }
#endif


#if HAVE_ALBERTA
    // albertagrid -> simplex
    template <class GridType = HostGrid,
      REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
        return MacroGridFactory<GridType>::createSimplexGrid(lower, upper, numCells);
      });
    }
#endif

    // yasp grid -> cube
    template <class GridType = HostGrid,
      class = typename GridType::YGrid>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
    {
      int overlap = 0;
      Parameters::get(name_ + "->overlap", overlap);
      std::string periodic(dimension, '0');
      Parameters::get(name_ + "->periodic", periodic); // e.g. 000 01 111

      return create_yaspgrid(Types<GridType>{}, overlap, std::bitset<dimension>(periodic));
    }

    template <int dim, class ct>
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>,
                                          int overlap, std::bitset<dimension> const& per) const
    {
      return create_structured_grid<int>([&](auto&& /*lower*/, auto&& upper, std::array<int,dimension> const& numCells)
      {
        return std::make_unique<HostGrid>(upper, numCells, per, overlap);
      });
    }

    template <int dim, class ct>
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>,
                                          int overlap, std::bitset<dimension> const& per) const
    {
      return create_structured_grid<int>([&](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
      {
        return std::make_unique<HostGrid>(lower, upper, numCells, per, overlap);
      });
    }

    template <int dim, class ct>
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>,
                                          int overlap, std::bitset<dimension> const& per) const
    {
      error_exit("MeshCreator cannot create YaspGrid with TensorProductCoordinates.");
      return {};
    }


#if HAVE_DUNE_SPGRID
    // spgrid -> cube
    template <class GridType = HostGrid,
      class = typename GridType::ReferenceCube,
      class = typename GridType::MultiIndex>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
    {
      return create_structured_grid<int>([](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
      {
        typename GridType::MultiIndex multiIndex(numCells);
        return std::make_unique<GridType>(lower, upper, multiIndex);
      });
    }
#endif

    // final fallback
    template <class GridType = HostGrid>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<0>) const
    {
      error_exit("Don't know how to create the grid.");
      return {};
    }

  private:
    std::string name_;
    mutable std::vector<int> boundaryIds_;
    mutable std::vector<int> elementIds_;
  };


} // end namespace AMDiS