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>
 
#if HAVE_DUNE_VTK
#include <dune/vtk/vtkreader.hh>
#include <dune/vtk/gridcreators/lagrangegridcreator.hh>
#endif
 
#if HAVE_DUNE_GMSH4
#include <dune/gmsh4/gmsh4reader.hh>
#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
#include <dune/gmsh4/utility/version.hh>
#endif
 
#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>
 
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
{
  template <class G>
  struct MeshCreator
  {
    enum { dimension = G::dimension };
    enum { dimworld = G::dimensionworld };
 
    using Grid = AdaptiveGrid_t<G>;
    using HostGrid = typename Grid::HostGrid;
 
    using ctype = typename Grid::ctype;
 
 
    MeshCreator(std::string const& name)
      : name_(name)
    {}
 
    static std::shared_ptr<Grid> create(std::string name)
    {
      return MeshCreator{name}.create();
    }
 
 
    std::unique_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.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));
    }
 
    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);
      });
    }
 
    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);
      });
    }
 
    std::vector<int> const& boundaryIds() const
    {
      return boundaryIds_;
    }
 
    std::vector<int> const& elementIds() const
    {
      return elementIds_;
    }
 
  private:
    static std::unique_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid)
    {
      return std::move(hostGrid);
    }
 
    template <class HG>
    static std::unique_ptr<Grid> construct(std::unique_ptr<HG> hostGrid)
    {
      return std::make_unique<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
    {
      if constexpr (int(Grid::dimension) == int(Grid::dimensionworld)) {
        // 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);
      } else {
        error_exit("Structured grids can be created for dim == dow only.");
        return nullptr;
      }
    }
 
    // 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") {
#if HAVE_DUNE_GMSH4
        if (Dune::Gmsh4::fileVersion(filename)[0] >= 4)
          return Dune::Gmsh4Reader<HostGrid, Dune::Gmsh4::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
        else
#endif
        return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{});
      }
#if HAVE_DUNE_VTK
      else if (ext == ".vtu") {
        return Dune::Vtk::VtkReader<HostGrid, Dune::Vtk::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
      }
#endif
      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> >>()) );
 
    template <class GridType, class LC>
    using SupportsInsertionIndex = decltype(std::declval<Dune::GridFactory<GridType>>().insertionIndex(std::declval<LC>()));
 
    // 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::GridFactory<GridType> factory;
      std::vector<int> boundaryIds, elementIds;
      Dune::GmshReader<GridType>::read(factory, filename, boundaryIds, elementIds);
 
      using HasInsertionIndexEntity
        = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::template Codim<0>::Entity>;
      using HasInsertionIndexIntersection
        = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::LeafIntersection>;
 
      auto gridPtr = factory.createGrid();
      if ((boundaryIds.empty() && elementIds.empty()) ||
          (!HasInsertionIndexEntity::value && !HasInsertionIndexIntersection::value))
        return std::unique_ptr<GridType>(std::move(gridPtr));
 
      // map boundaryIds and elementIds read from file to grid indexing.
      if (!boundaryIds.empty() && HasInsertionIndexIntersection::value)
        boundaryIds_.resize(gridPtr->numBoundarySegments());
      if (!elementIds.empty() && HasInsertionIndexEntity::value)
        elementIds_.resize(gridPtr->size(0));
 
      auto const& indexSet = gridPtr->leafIndexSet();
      for (auto const& e : elements(gridPtr->leafGridView())) {
        if constexpr(HasInsertionIndexEntity::value) {
          if (!elementIds.empty())
            elementIds_[indexSet.index(e)] = elementIds[factory.insertionIndex(e)];
        }
 
        if (boundaryIds.empty() || !e.hasBoundaryIntersections())
          continue;
 
        for (auto const& it : intersections(gridPtr->leafGridView(), e)) {
          if constexpr(HasInsertionIndexIntersection::value) {
            if (it.boundary())
              boundaryIds_[it.boundarySegmentIndex()]
                = boundaryIds[factory.insertionIndex(it)];
          }
        }
      }
 
      return std::unique_ptr<GridType>(std::move(gridPtr));
    }
 
    // fallback if GmshReader cannot be used
    template <class GridType = HostGrid>
    std::unique_ptr<GridType> read_gmsh_file(std::string const&, 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&, Dune::PriorityTag<0>) const
    {
      error_exit("AlbertaGrid (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
      return {};
    }
#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, std::bitset<dimension> const&) 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