From 783a01f17b1fdf8fa8251e6572dbf2d1a2c82476 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Sat, 1 Jun 2019 15:47:56 +0200
Subject: [PATCH] [subgrid] Restructure. Add option for reading from binary
file
Based on initial work by Kilian Weishaupt.
---
dumux/io/grid/subgridmanager.hh | 266 ++++++++++++++++++++++++++++----
1 file changed, 236 insertions(+), 30 deletions(-)
diff --git a/dumux/io/grid/subgridmanager.hh b/dumux/io/grid/subgridmanager.hh
index e5cd7e258e..01523b8b9a 100644
--- a/dumux/io/grid/subgridmanager.hh
+++ b/dumux/io/grid/subgridmanager.hh
@@ -28,66 +28,101 @@
#include <memory>
+#include <dune/common/deprecated.hh>
+#include <dune/common/shared_ptr.hh>
+#include <dune/common/concept.hh>
#include <dune/subgrid/subgrid.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/boundaryflag.hh>
+#include <dumux/io/grid/gridmanager_base.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
-// TODO: remove this after 3.1 is released
+#include <dumux/io/rasterimagereader.hh>
+
+// TODO: deprecated: remove this after 3.1 is released
#include <dune/grid/io/file/vtk.hh>
-// TODO: remove this after 3.1 is released
+// TODO: deprecated: remove this after 3.1 is released
#include <dune/grid/io/file/dgfparser/dgfwriter.hh>
namespace Dumux {
+namespace Concept {
+
+/*!
+ * \ingroup InputOutput
+ * \brief The element selector concept
+ */
+template<class Element>
+struct ElementSelector
+{
+ template<class F>
+ auto require(F&& f) -> decltype(
+ bool(f(std::declval<const Element&>()))
+ );
+};
+} // end namespace Concept
/*!
* \ingroup InputOutput
- * \brief A grid manager for dune-subgrid.
+ * \brief The base class for grid managers for dune-subgrid.
*/
-template <class HostGrid>
-class SubgridManager
+template <class HostGrid, class HostGridManager>
+class SubGridManagerBase
{
- static constexpr auto dim = HostGrid::dimension;
+ static constexpr int dim = HostGrid::dimension;
+ using HostElement = typename HostGrid::template Codim<0>::Entity;
public:
using Grid = Dune::SubGrid<dim, HostGrid>;
+ /*!
+ * \brief Make the grid using an externally created host grid.
+ */
+ template<class ES,
+ typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0>
+ void init(HostGrid& hostGrid,
+ const ES& selector,
+ const std::string& paramGroup = "")
+ {
+ subGrid_ = createGrid_(hostGrid, selector, paramGroup);
+ loadBalance();
+ }
+
+ /*!
+ * \brief Make the grid and create the host grid internally.
+ */
+ template<class ES,
+ typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0>
+ void init(const ES& selector,
+ const std::string& paramGroup = "")
+ {
+ initHostGrid_(paramGroup);
+ subGrid_ = createGrid_(hostGridManager_->grid(), selector, paramGroup);
+ loadBalance();
+ }
+
/*!
* \brief Make the subgrid.
*/
template<class ElementSelector>
- static std::unique_ptr<Grid> makeGrid(HostGrid& hostgrid,
+ DUNE_DEPRECATED_MSG("Create an instance of this class and use subgridManager.init(hostGrid, selector, paramGroup)")
+ static std::unique_ptr<Grid> makeGrid(HostGrid& hostGrid,
const ElementSelector& selector,
- const std::string& modelParamGroup = "")
+ const std::string& paramGroup = "")
{
- // A unique pointer to the subgrid.
- auto subgridPtr = std::make_unique<Grid>(hostgrid);
+ std::cerr << "Deprecation warning: SubGridManager::makeGrid is deprecated."
+ << "Create an instance of this class and use subgridManager.init(hostGrid, selector, paramGroup)." << std::endl;
- // A container to store the host grid elements' ids.
- std::set<typename HostGrid::Traits::GlobalIdSet::IdType> elementsForSubgrid;
- const auto& globalIDset = subgridPtr->getHostGrid().globalIdSet();
-
- // Construct the subgrid.
- subgridPtr->createBegin();
-
- // Loop over all elements of the host grid and use the selector to
- // choose which elements to add to the subgrid.
- auto hostGridView = subgridPtr->getHostGrid().leafGridView();
- for (const auto& e : elements(hostGridView))
- if (selector(e))
- elementsForSubgrid.insert(globalIDset.template id<0>(e));
-
- subgridPtr->insertSetPartial(elementsForSubgrid);
- subgridPtr->createEnd();
+ auto subgridPtr = createGrid_(hostGrid, selector, paramGroup);
// TODO: remove this after 3.1 is released
// If desired, write out the final subgrid as a dgf file.
- if (getParamFromGroup<bool>(modelParamGroup, "Grid.WriteSubGridToDGF", false))
+ if (getParamFromGroup<bool>(paramGroup, "Grid.WriteSubGridToDGF", false))
{
std::cerr << "Deprecation warning: SubGridManager: Grid.WriteSubGridToDGF is deprecated."
<< "Use Dune::VTKWriter to write out your grid manually." << std::endl;
- const auto postfix = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name", "");
+ const auto postfix = getParamFromGroup<std::string>(paramGroup, "Problem.Name", "");
const std::string name = postfix == "" ? "subgrid" : "subgrid_" + postfix;
Dune::DGFWriter<typename Grid::LeafGridView> writer(subgridPtr->leafGridView());
writer.write(name + ".dgf");
@@ -95,12 +130,12 @@ public:
// TODO: remove this after 3.1 is released
// If desired, write out the hostgrid as vtk file.
- if (getParamFromGroup<bool>(modelParamGroup, "Grid.WriteSubGridToVtk", false))
+ if (getParamFromGroup<bool>(paramGroup, "Grid.WriteSubGridToVtk", false))
{
std::cerr << "Deprecation warning: SubGridManager: Grid.WriteSubGridToVtk is deprecated."
<< "Use Dune::VTKWriter to write out your grid manually." << std::endl;
- const auto postfix = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name", "");
+ const auto postfix = getParamFromGroup<std::string>(paramGroup, "Problem.Name", "");
const std::string name = postfix == "" ? "subgrid" : "subgrid_" + postfix;
Dune::VTKWriter<typename Grid::LeafGridView> vtkWriter(subgridPtr->leafGridView());
vtkWriter.write(name);
@@ -109,6 +144,177 @@ public:
// Return a unique pointer to the subgrid.
return subgridPtr;
}
+
+ /*!
+ * \brief Returns a reference to the grid.
+ */
+ Grid& grid()
+ {
+ return *subGrid_;
+ }
+
+ /*!
+ * \brief Distributes the grid on all processes of a parallel
+ * computation.
+ */
+ void loadBalance()
+ {
+ subGrid_->loadBalance();
+ }
+
+protected:
+
+ /*!
+ * \brief Make the subgrid.
+ */
+ template<class ES,
+ typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0>
+ static std::unique_ptr<Grid> createGrid_(HostGrid& hostGrid,
+ const ES& selector,
+ const std::string& paramGroup = "")
+ {
+ // A unique pointer to the subgrid.
+ auto subgridPtr = std::make_unique<Grid>(hostGrid);
+
+ // A container to store the host grid elements' ids.
+ std::set<typename HostGrid::Traits::GlobalIdSet::IdType> elementsForSubgrid;
+ const auto& globalIDset = subgridPtr->getHostGrid().globalIdSet();
+
+ // Construct the subgrid.
+ subgridPtr->createBegin();
+
+ // Loop over all elements of the host grid and use the selector to
+ // choose which elements to add to the subgrid.
+ auto hostGridView = subgridPtr->getHostGrid().leafGridView();
+ for (const auto& e : elements(hostGridView))
+ if (selector(e))
+ elementsForSubgrid.insert(globalIDset.template id<0>(e));
+
+ if (elementsForSubgrid.empty())
+ DUNE_THROW(Dune::GridError, "No elements in subgrid");
+
+ subgridPtr->insertSetPartial(elementsForSubgrid);
+ subgridPtr->createEnd();
+
+ // Return a unique pointer to the subgrid.
+ return subgridPtr;
+ }
+
+ void initHostGrid_(const std::string& paramGroup)
+ {
+ hostGridManager_ = std::make_unique<HostGridManager>();
+ hostGridManager_->init(paramGroup);
+ }
+
+ /*!
+ * \brief Returns a reference to the host grid.
+ */
+ HostGrid& hostGrid_()
+ {
+ return hostGridManager_->grid();
+ }
+
+ std::shared_ptr<Grid> subGrid_;
+ std::unique_ptr<HostGridManager> hostGridManager_;
+};
+
+/*!
+ * \ingroup InputOutput
+ * \brief A grid manager for dune-subgrid
+ */
+template<class HostGrid, class HostGridManager = GridManager<HostGrid>>
+class SubgridManager : public SubGridManagerBase<HostGrid, HostGridManager>
+{};
+
+/*!
+ * \ingroup InputOutput
+ * \brief A grid manager for dune-subgrid using YaspGrid as the host grid
+ * With YaspGrid as host grid we also support creating grids from raster images
+ */
+template<int dim, class Coordinates>
+class SubgridManager<Dune::YaspGrid<dim, Coordinates>>
+: public SubGridManagerBase<Dune::YaspGrid<dim, Coordinates>, GridManager<Dune::YaspGrid<dim, Coordinates>>>
+{
+ using ParentType = SubGridManagerBase<Dune::YaspGrid<dim, Coordinates>, GridManager<Dune::YaspGrid<dim, Coordinates>>>;
+public:
+ using ParentType::Grid;
+ using ParentType::init;
+
+ /*!
+ * \brief Make the subgrid without host grid and element selector
+ * This means we try to construct the element selector from the input file
+ * \param paramGroup the parameter file group to check
+ */
+ void init(const std::string& paramGroup = "")
+ {
+ // check if there is an image file we can construct the element selector from
+ if (hasParamInGroup(paramGroup, "Grid.Image"))
+ {
+ const auto imgFileName = getParamFromGroup<std::string>(paramGroup, "Grid.Image");
+ const auto ext = getFileExtension_(imgFileName);
+ if (ext == "pbm")
+ {
+ if (dim != 2)
+ DUNE_THROW(Dune::GridError, "Portable Bitmap Format only supports dim == 2");
+
+ // read image
+ const auto img = NetPBMReader::readPBM(imgFileName);
+ createGridFromImage_(img, paramGroup);
+ }
+ else
+ DUNE_THROW(Dune::IOError, "The SubGridManager doesn't support image files with extension: *." << ext);
+
+ }
+ else
+ DUNE_THROW(Dune::IOError, "SubGridManager couldn't construct element selector. Specify Grid.Image in the input file!");
+ }
+
+private:
+ template<class Img>
+ void createGridFromImage_(const Img& img, const std::string& paramGroup)
+ {
+ // check if the number of cells matches
+ std::array<int, dim> cells; cells.fill(1);
+ cells = getParamFromGroup<std::array<int, dim>>(paramGroup, "Grid.Cells", cells);
+
+ if (img.header().nCols != cells[0])
+ DUNE_THROW(Dune::GridError, "Host grid has wrong number of cells in x-direction. Expected "
+ << img.header().nCols << ", got " << cells[0]);
+ if (img.header().nRows != cells[1])
+ DUNE_THROW(Dune::GridError, "Host grid has wrong number of cells in y-direction. Expected "
+ << img.header().nRows << ", got " << cells[1]);
+
+
+ // construct the host grid
+ this->initHostGrid_(paramGroup);
+
+ // check if the marker is customized, per default
+ // we mark all cells that are encoded as 0
+ const bool marked = getParamFromGroup<bool>(paramGroup, "Grid.Marker", 0);
+
+ // Create the selector
+ auto elementSelector = [this, &img, &marked](const auto& element)
+ {
+ const auto eIdx = this->hostGrid_().leafGridView().indexSet().index(element);
+ return img[eIdx] == marked;
+ };
+
+ // create the grid
+ this->subGrid_ = this->createGrid_(this->hostGrid_(), elementSelector, paramGroup);
+ this->loadBalance();
+ }
+
+ /*!
+ * \brief Returns the filename extension of a given filename
+ */
+ std::string getFileExtension_(const std::string& fileName) const
+ {
+ std::size_t i = fileName.rfind('.', fileName.length());
+ if (i != std::string::npos)
+ return(fileName.substr(i+1, fileName.length() - i));
+ else
+ DUNE_THROW(Dune::IOError, "Please provide and extension for your grid file ('"<< fileName << "')!");
+ }
};
//! dune-subgrid doesn't have this implemented
--
GitLab