Mesh.hpp

#pragma once

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

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

#include <dune/grid/albertagrid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>

#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/gmshreader.hh>

#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/Output.hpp>
#include <dune/amdis/utility/Filesystem.hpp>

namespace AMDiS
{
  namespace tag
  {
    struct albertagrid {};
    struct uggrid {};
    struct yaspgrid {};
    struct unknowngrid {};

  } // end namespace tag

  namespace Impl
  {
    template <class Grid>
    struct MeshTag
    {
      using type = tag::unknowngrid;
    };

    // specialization for some grid types from DUNE
#if HAVE_ALBERTA
    template <int dim, int dimworld>
    struct MeshTag<Dune::AlbertaGrid<dim, dimworld>>
    {
      using type = tag::albertagrid;
    };
#endif

#if HAVE_UG
    template <int dim>
    struct MeshTag<Dune::UGGrid<dim>>
    {
      using type = tag::uggrid;
    };
#endif

    template <int dim, class Coordinates>
    struct MeshTag<Dune::YaspGrid<dim, Coordinates>>
    {
      using type = tag::yaspgrid;
    };

  } // end namespace Impl

  template <class Grid>
  using MeshTag_t = typename Impl::MeshTag<Grid>::type;


  /// A creator class for meshes. Each mesh needs different way of initialization
  template <class Grid>
  class MeshCreator
  {
    static std::unique_ptr<Grid> create(std::string meshName)
    {
      error_exit("Creator not yet implemented for this mesh type.");
    }
  };

#if HAVE_ALBERTA
  template <int dim, int dimworld>
  struct MeshCreator<Dune::AlbertaGrid<dim, dimworld>>
  {
    using Grid = Dune::AlbertaGrid<dim, dimworld>;

    static std::unique_ptr<Grid> create(std::string meshName)
    {
      std::string macro_filename = "";
      Parameters::get(meshName + "->macro file name", macro_filename);

      // TODO: if filename_extension is ".2d" or ".3d" read it directly from file
      // otherwise use a factory method

      return std::make_unique<Grid>(macro_filename);
    }
  };
#endif


#if HAVE_UG
  template <int dim>
  struct MeshCreator<Dune::UGGrid<dim>>
  {
    using Grid = Dune::UGGrid<dim>;

    static std::unique_ptr<Grid> create(std::string meshName)
    {

      std::string filename = "";
      Parameters::get(meshName + "->macro file name", filename);

      if (!filename.empty()) {
        filesystem::path fn(filename);
        auto ext = fn.extension();

#if HAVE_ALBERTA
        if (ext == "1d" || ext == "2d" || ext == "3d") {
          Dune::GridFactory<Grid> factory;
          Dune::AlbertaReader<Grid> reader;
          reader.readGrid(filename, factory);
          return std::unique_ptr<Grid>{factory.createGrid()};
        }
#endif
        if (ext == "msh") {
          Dune::GmshReader<Grid> reader;
          return std::unique_ptr<Grid>{reader.read(filename)};
        }
      } else {
        error_exit("Construction of UGGrid without filename not yet implemented!");
      }

      error_exit("No way to construct UG-Grid found");
      return {};
    }
  };
#endif

  template <int dim, class T>
  struct MeshCreator<Dune::YaspGrid<dim, Dune::EquidistantCoordinates<T,dim>>>
  {
    using Grid = Dune::YaspGrid<dim, Dune::EquidistantCoordinates<T,dim>>;

    static std::unique_ptr<Grid> create(std::string meshName)
    {
      Dune::FieldVector<double, dim> L; L = 1.0;  // extension of the domain
      Parameters::get(meshName + "->dimension", L);

      msg("L = ", L);

      auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
      Parameters::get(meshName + "->num cells", s);

      // TODO: add more parameters for yasp-grid (see constructor)

      return std::make_unique<Grid>(L, s);
    }
  };


  template <int dim, class T>
  struct MeshCreator<Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<T, dim>>>
  {
    using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<T, dim>>;

    static std::unique_ptr<Grid> create(std::string meshName)
    {
      Dune::FieldVector<double, dim> lowerleft;  lowerleft = 0.0; // Lower left corner of the domain
      Dune::FieldVector<double, dim> upperright; upperright = 1.0; // Upper right corner of the domain
      Parameters::get(meshName + "->min corner", lowerleft);
      Parameters::get(meshName + "->max corner", upperright);

      auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
      Parameters::get(meshName + "->num cells", s);

      // TODO: add more parameters for yasp-grid (see constructor)

      return std::make_unique<Grid>(lowerleft, upperright, s);
    }
  };

} // end namespace AMDiS