rewritten mesh creator to support more general creation

examples/cahn_hilliard.cc

@@ -6,6 +6,8 @@
 #include <amdis/GridFunctions.hpp>
 #include <amdis/Marker.hpp>
 
+#include <dune/grid/albertagrid.hh>
+
 using namespace AMDiS;
 
 using Grid = Dune::AlbertaGrid<GRIDDIM, WORLDDIM>;

examples/init/heat.dat.2d

 heatMesh->macro file name:    ./macro/macro.stand.2d
 heatMesh->global refinements: 4
 heatMesh->min corner: 0 0
-heatMesh->max corner: 1 1
+heatMesh->max corner: 2 2
 heatMesh->num cells: 2 2
-heatMesh->dimension: 2 2
 
 heat->mesh:  heatMesh
 heat->names: u

examples/init/stokes.dat.2d

 stokesMesh->global refinements: 0
-stokesMesh->dimension: 1.0 1.0
+stokesMesh->max corner: 1.0 1.0
 stokesMesh->num cells: 4 4
 
 stokes->mesh: stokesMesh

src/amdis/BoundaryManager.hpp

@@ -154,6 +154,13 @@ namespace AMDiS
       }
     }
 
+    template <class I>
+    void setBoundaryIds(std::vector<I> const& ids)
+    {
+      test_exit(ids.size() == boundaryIds_.size(), "Number of boundary IDs does not match!");
+      std::copy(ids.begin(), ids.end(), boundaryIds_.begin());
+    }
+
   private:
     std::shared_ptr<Grid> grid_;
     using Super::boundaryIds_;

src/amdis/CMakeLists.txt

@@ -50,7 +50,7 @@ install(FILES
     LocalOperators.hpp
     Marker.hpp
     Marker.inc.hpp
-    Mesh.hpp
+    MeshCreator.hpp
     Operations.hpp
     OperatorList.hpp
     Output.hpp

src/amdis/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 <amdis/Initfile.hpp>
-#include <amdis/Output.hpp>
-#include <amdis/common/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>
-  struct MeshCreator
-  {
-    static std::unique_ptr<Grid> create(std::string const& meshName)
-    {
-      error_exit("Creator not yet implemented for this mesh type.");
-      return {};
-    }
-  };
-
-#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 const& 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 const& meshName)
-    {
-
-      std::string filename = "";
-      Parameters::get(meshName + "->macro file name", filename);
-
-      test_exit(!filename.empty(),
-        "Construction of UGGrid without filename not yet implemented!");
-
-      filesystem::path fn(filename);
-      auto ext = fn.extension();
-
-      if (ext == ".msh") {
-        Dune::GmshReader<Grid> reader;
-        return std::unique_ptr<Grid>{reader.read(filename)};
-      }
-// #if HAVE_ALBERTA
-//       else if (ext == ".1d" || ext == ".2d" || ext == ".3d") {
-//         Dune::Hybrid::ifElse(bool_t<ALBERTA_DIM == Grid::dimensionworld>{},
-//         [&](auto id) {
-//           Dune::GridFactory<Grid> factory;
-//           Dune::AlbertaReader<Grid> reader;
-//           id(reader).readGrid(filename, id(factory));
-//           return std::unique_ptr<Grid>{factory.createGrid()};
-//         });
-
-//         warning("WORLDDIM must be given in Alberta flags.");
-//       }
-// #endif
-
-      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 const& meshName)
-    {
-      Dune::FieldVector<T, dim> L; L = 1.0;  // extension of the domain
-      Parameters::get(meshName + "->dimension", 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 const& meshName)
-    {
-      Dune::FieldVector<T, dim> lowerleft;  lowerleft = 0.0; // Lower left corner of the domain
-      Dune::FieldVector<T, 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

src/amdis/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/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 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;
+
+    /// Construct a new MeshCreator
+    /**
+     * \param name  The name of the mesh used in the initifile
+     **/
+    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
+     *
+     * Otherwise tries to create a grid depending on the grid type.
+     **/
+    std::unique_ptr<Grid> create() const
+    {
+      auto filename = Parameters::get<std::string>(name_ + "->macro file name");
+      auto structured = Parameters::get<std::string>(name_ + "->structured");
+      if (filename) {
+        // read a macro file
+        return create_unstructured_grid(filename.value());
+      } else if (structured) {
+        // create structured grids
+        if (structured.value() == "cube") {
+          return create_cube_grid();
+        } else if (structured && structured.value() == "simplex") {
+          return create_simplex_grid();
+        } else {
+          error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
+          return {};
+        }
+      } else {
+        // decide by inspecting the grid type how to create the grid
+        return create_by_gridtype<Grid>(Dune::PriorityTag<42>{});
+      }
+    }
+
+    /// Create a structured cube grid
+    std::unique_ptr<Grid> create_cube_grid() const
+    {
+      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
+      {
+        return Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, numCells);
+      });
+    }
+
+    /// Create a structured simplex grid
+    std::unique_ptr<Grid> create_simplex_grid() const
+    {
+      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
+      {
+        return MacroGridFactory<Grid>::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:
+    // use the structured grid factory to create the grid
+    template <class Factory>
+    std::unique_ptr<Grid> 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),unsigned int>(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<Grid> create_unstructured_grid(std::string const& filename) const
+    {
+      filesystem::path fn(filename);
+      auto ext = fn.extension();
+
+      if (ext == ".msh") {
+        return read_gmsh_file<Grid>(filename, Dune::PriorityTag<42>{});
+      }
+      else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
+        return read_alberta_file<Grid>(filename, Dune::PriorityTag<42>{});
+      }
+      else {
+        error_exit("Can not 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 = Grid,
+      REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)>
+    std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<2>) const
+    {
+      Dune::GmshReader<GridType> reader;
+      return std::unique_ptr<GridType>{reader.read(filename, boundaryIds_, elementIds_)};
+    }
+
+    // fallback if GmshReader cannot be used
+    template <class GridType = Grid>
+    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 = Grid,
+      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 = Grid,
+      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;
+      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 = Grid,
+      REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
+    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
+    {
+      error_exit("Alberta (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("Alberta (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
+      return {};
+    }
+#endif
+
+
+#if HAVE_ALBERTA
+    // albertagrid -> simplex
+    template <class GridType = Grid,
+      REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
+    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const
+    {
+      return create_simplex_grid();
+    }
+#endif
+
+    // yasp grid -> cube
+    template <class GridType = Grid,
+      class = typename GridType::YGrid>
+    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
+    {
+      return create_cube_grid();
+    }
+
+#if HAVE_DUNE_SPGRID
+    // spgrid -> cube
+    template <class GridType = Grid,
+      class = typename GridType::ReferenceCube,
+      class = typename GridType::MultiIndex>
+    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
+    {
+      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
+      {
+        std::array<int,dimworld> cells;
+        std::copy(std::begin(numCells), std::end(numCells), std::begin(cells));
+        typename GridType::MultiIndex multiIndex(cells);
+        return std::make_unique<GridType>(lower, upper, multiIndex);
+      });
+    }
+#endif
+
+    // final fallback
+    template <class GridType = Grid>
+    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

src/amdis/ProblemStat.hpp

@@ -22,7 +22,7 @@
 #include <amdis/LinearAlgebra.hpp>
 #include <amdis/OperatorList.hpp>
 #include <amdis/Marker.hpp>
-#include <amdis/Mesh.hpp>
+#include <amdis/MeshCreator.hpp>
 #include <amdis/PeriodicBC.hpp>
 #include <amdis/ProblemStatBase.hpp>
 #include <amdis/ProblemStatTraits.hpp>

src/amdis/ProblemStat.inc.hpp

@@ -126,8 +126,11 @@ template <class Traits>
 void ProblemStat<Traits>::createGrid()
 {
   Parameters::get(name_ + "->mesh", gridName_);
-  grid_ = MeshCreator<Grid>::create(gridName_);
+  MeshCreator<Grid> creator(gridName_);
+  grid_ = creator.create();
   boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
+  if (!creator.boundaryIds().empty())
+    boundaryManager_->setBoundaryIds(creator.boundaryIds());
 
   msg("Create grid:");
   msg("#elements = {}"   , grid_->size(0));

src/amdis/common/Concepts.hpp

@@ -53,6 +53,17 @@ namespace AMDiS
     struct IsReferenceWrapper<std::reference_wrapper<T>>
       : std::true_type {};
 
+
+    template <class, class = void>
+    struct IsDefined
+      : std::false_type {};
+
+    template <class T>
+    struct IsDefined<T, std::enable_if_t<std::is_object<T>::value &&
+                                        !std::is_pointer<T>::value &&
+                                        (sizeof(T) > 0)> >
+      : std::true_type {};
+
   } // end namespace Traits
 
   namespace Concepts
@@ -117,6 +128,9 @@ namespace AMDiS
     template <class MI>
     constexpr bool MultiIndex = models<Definition::MultiIndex(MI)>;
 
+    template <class T>
+    constexpr bool Defined = Traits::IsDefined<T>::value;
+
     /** @} **/
 
   } // end namespace Concepts

src/amdis/utility/CMakeLists.txt

@@ -2,5 +2,6 @@ install(FILES
     AssembleOperators.hpp
     LocalBasisCache.hpp
     LocalToGlobalAdapter.hpp
+    MacroGridFactory.hpp
     QuadratureFactory.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility)

src/amdis/utility/MacroGridFactory.hpp

+#pragma once
+
+#include <dune/common/typeutilities.hh>
+#include <dune/common/version.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#if ! DUNE_VERSION_GT(DUNE_GRID,2,6)
+  #include <dune/grid/yaspgrid.hh>
+#endif
+
+namespace Dune
+{
+  template <class GridType>
+  struct CubedWrapper : public GridType {};
+
+
+  template <class GridType>
+  class GridFactory<CubedWrapper<GridType>>
+      : public GridFactoryInterface<GridType>
+  {
+    using Self = GridFactory;
+    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))
+    {}
+
+    /// \brief Insert a vertex into the coarse grid
+    void insertVertex (const FieldVector<ctype,dimworld>& pos) override
+    {
+      factory_->insertVertex(pos);
+    }
+
+    /// \brief Insert simplex elements into the coarse grid
+    /**
+     * Creates a simplex subdividion of the cube element.
+     *
+     * \param type      The GeometryType of the box grid
+     * \param vertices  Indices of the cube corners
+     **/
+    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