Merge branch 'feature/implement-gmshreader' into 'master'

Enable GmshReader for multiple domains

See merge request !43

dune.module

@@ -7,4 +7,5 @@ Version: 2.8-git
 Maintainer: santiago.ospina@iwr.uni-heidelberg.de
 #depending on
 Depends: dune-grid (>= 2.8)
+Suggests: dune-uggrid
 Whitespace-Hook: Yes

dune/grid/multidomaingrid.hh

@@ -4,6 +4,8 @@
 #include <dune/common/parallel/communication.hh>
 #include <dune/grid/multidomaingrid/multidomaingrid.hh>
 #include <dune/grid/multidomaingrid/multidomainmcmgmapper.hh>
+#include <dune/grid/multidomaingrid/factory.hh>
+#include <dune/grid/multidomaingrid/gmshreader.hh>
 
 namespace Dune {
 

dune/grid/multidomaingrid/factory.hh

+#ifndef DUNE_MULTIDOMAINGRID_FACTORY_HH
+#define DUNE_MULTIDOMAINGRID_FACTORY_HH
+
+namespace Dune {
+
+namespace mdgrid {
+
+// forward declaration
+template<typename HostGrid,typename MDGridTraits>
+class MultiDomainGrid;
+
+} // namespace mdgrid
+
+// forward declaration
+
+template <class GridType>
+class GridFactoryInterface;
+
+template<class HostGrid, class MDGTraits>
+class GridFactory<Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>> : public GridFactoryInterface<Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>> {
+
+  using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+  using Base = GridFactoryInterface<Grid>;
+public:
+
+  GridFactory()
+    : _host_grid_factory{}
+  {}
+
+  // use default implementation from base class
+  using Base::insertBoundarySegment;
+
+  void insertVertex(const FieldVector<typename Grid::ctype, Grid::dimensionworld>& pos) override {
+    _host_grid_factory.insertVertex(pos);
+  }
+
+  void insertElement(const GeometryType& type, const std::vector<unsigned int>& vertices) override {
+    _host_grid_factory.insertElement(type, vertices);
+  }
+
+  void insertBoundarySegment(const std::vector<unsigned int>& vertices) override {
+    _host_grid_factory.insertBoundarySegment(vertices);
+  }
+
+  std::unique_ptr<Grid> createGrid() override  {
+    if (not _grid_ptr)
+      makeGrid();
+    assert(_grid_ptr);
+    return std::exchange(_grid_ptr, nullptr);
+  }
+
+  void makeGrid(int max_subdomains = 1) {
+    assert(not _grid_ptr);
+
+    std::unique_ptr<MDGTraits> traits;
+    if constexpr (std::is_default_constructible_v<MDGTraits>)
+      traits = std::make_unique<MDGTraits>();
+    else
+      traits = std::make_unique<MDGTraits>(max_subdomains);
+
+    if (traits->maxSubDomainIndex() < max_subdomains)
+      DUNE_THROW(IOError, "Maximum number of sub-domains is bigger than what grid traits allows");
+
+    _grid_ptr = std::make_unique<Grid>(_host_grid_factory.createGrid(), *traits);
+  }
+
+  Grid& grid() {
+    assert(_grid_ptr);
+    return *_grid_ptr;
+  }
+
+  HostGrid& hostGrid() {
+    assert(_host_grid_ptr);
+    return *_host_grid_ptr;
+  }
+
+  const std::shared_ptr<HostGrid>& hostGridPtr() {
+    return _host_grid_ptr;
+  }
+
+
+  Dune::GridFactory<HostGrid>& hostGridFactory() {
+    return _host_grid_factory;
+  }
+
+private:
+  Dune::GridFactory<HostGrid> _host_grid_factory;
+  std::shared_ptr<HostGrid> _host_grid_ptr;
+  std::unique_ptr<Grid> _grid_ptr;
+};
+
+} // namespace Dune
+
+#endif // DUNE_MULTIDOMAINGRID_FACTORY_HH

dune/grid/multidomaingrid/gmshreader.hh

+#ifndef DUNE_MULTIDOMAINGRID_GMSHREADERPARSER_HH
+#define DUNE_MULTIDOMAINGRID_GMSHREADERPARSER_HH
+
+namespace Dune {
+
+namespace mdgrid {
+
+// forward declaration
+template<typename HostGrid,typename MDGridTraits>
+class MultiDomainGrid;
+
+} // namespace mdgrid
+
+// forward declaration
+template<typename GridType>
+class GmshReaderParser;
+
+template<class HostGrid, class MDGTraits>
+class GmshReaderParser<Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>> : public Dune::GmshReaderParser<HostGrid> {
+public:
+
+  using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+
+  GmshReaderParser(Dune::GridFactory<Grid>& factory, bool v, bool i)
+    : Dune::GmshReaderParser<HostGrid>{factory.hostGridFactory(), v, i}
+    , _factory{factory}
+  {}
+
+  void read (const std::string& f)
+  {
+    Dune::GmshReaderParser<HostGrid>::read(f);
+    const auto& index_map = this->elementIndexMap();
+    int max_subdomains = *std::max_element(begin(index_map), end(index_map));
+    _factory.makeGrid(max_subdomains);
+    auto& grid = _factory.grid();
+
+    grid.startSubDomainMarking();
+    unsigned int i = 0;
+    for (const auto& cell : elements(grid.leafGridView())) {
+      auto subdomain = index_map[i] - 1; // gmsh index starts at 1!
+      grid.addToSubDomain(subdomain, cell), i++;
+    }
+
+    grid.preUpdateSubDomains();
+    grid.updateSubDomains();
+    grid.postUpdateSubDomains();
+  }
+
+private:
+  Dune::GridFactory<Grid>& _factory;
+};
+
+} // namespace Dune
+
+#endif // DUNE_MULTIDOMAINGRID_GMSHREADERPARSER_HH

test/CMakeLists.txt

@@ -58,3 +58,9 @@ dune_add_test(
   )
 
 dune_add_test(SOURCES testpartitioning.cc)
+
+dune_add_test(
+  SOURCES testgmshreader.cc
+  CMD_ARGS "${CMAKE_CURRENT_SOURCE_DIR}/example.msh" 3
+  CMAKE_GUARD dune-uggrid_FOUND
+)

test/testgmshreader.cc

+
+#include "config.h"
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/multidomaingrid.hh>
+#include <iostream>
+
+int main(int argc, char** argv)
+{
+
+  try {
+
+    Dune::MPIHelper::instance(argc,argv);
+
+    if (argc != 3) {
+      std::cerr << "Usage: " << argv[0] << " <gmhs-file> <expected-domains>" << std::endl;
+      DUNE_THROW(Dune::IOError, "");
+    }
+
+    std::string fileName = argv[1];
+    std::size_t expectedDomains = atoi(argv[2]);
+
+    using HostGrid = Dune::UGGrid<2>;
+    using MDGrid = Dune::MultiDomainGrid<HostGrid,Dune::mdgrid::FewSubDomainsTraits<2,8> >;
+
+    auto mdgrid = Dune::GmshReader<MDGrid>::read(fileName);
+
+    if (mdgrid->maxAssignedSubDomainIndex() + 1 != expectedDomains)
+      DUNE_THROW(Dune::InvalidStateException, "Grid does not contain expected number of sub-domains");
+
+  } catch (Dune::Exception &e) {
+    std::cerr << e << std::endl;
+    return 1;
+  } catch (...) {
+    std::cerr << "Generic exception!" << std::endl;
+    return 2;
+  }
+
+}