diff --git a/dumux/io/grid/gridmanager_sub.hh b/dumux/io/grid/gridmanager_sub.hh
index 16f884e86fdbd2e3a3f0f4a4466e9210db45cc7e..cf18c0bb848ffe93732876667b8a01fe69435ed7 100644
--- a/dumux/io/grid/gridmanager_sub.hh
+++ b/dumux/io/grid/gridmanager_sub.hh
@@ -70,6 +70,7 @@ class SubGridManagerBase
{
static constexpr int dim = HostGrid::dimension;
using HostElement = typename HostGrid::template Codim<0>::Entity;
+ using GlobalPosition = typename HostElement::Geometry::GlobalCoordinate;
public:
using Grid = Dune::SubGrid<dim, HostGrid>;
@@ -153,6 +154,16 @@ protected:
hostGridManager_->init(paramGroup);
}
+ void initHostGrid_(const GlobalPosition& lowerLeft,
+ const GlobalPosition& upperRight,
+ const std::array<int, dim>& cells,
+ const std::string& paramGroup,
+ const int overlap = 1)
+ {
+ hostGridManager_ = std::make_unique<HostGridManager>();
+ hostGridManager_->init(lowerLeft, upperRight, cells, paramGroup, overlap);
+ }
+
/*!
* \brief Returns a reference to the host grid.
*/
@@ -230,20 +241,16 @@ 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]);
+ // get the corner coordinates
+ using GlobalPosition = typename ParentType::Grid::template Codim<0>::Geometry::GlobalCoordinate;
+ const auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.UpperRight");
+ const auto lowerLeft = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.LowerLeft", GlobalPosition(0.0));
+ // get the number of cells
+ const std::array<int, dim> cells{static_cast<int>(img.header().nCols), static_cast<int>(img.header().nRows)};
// construct the host grid
- this->initHostGrid_(paramGroup);
+ this->initHostGrid_(lowerLeft, upperRight, cells, paramGroup);
// check if the marker is customized, per default
// we mark all cells that are encoded as 0
@@ -252,12 +259,26 @@ private:
// Create the selector
auto elementSelector = [this, &img, &marked](const auto& element)
{
- const auto eIdx = this->hostGrid_().leafGridView().indexSet().index(element);
+ auto eIdx = this->hostGrid_().leafGridView().indexSet().index(element);
+
+ // if the hostgrid was refined, get the index of the original, un-refined
+ // host grid element which fits with the image's indices
+ if (element.hasFather())
+ {
+ auto level0Element = element.father();
+ while (level0Element.hasFather())
+ level0Element = level0Element.father();
+
+ assert(level0Element.level() == 0);
+ eIdx = this->hostGrid_().levelGridView(0).indexSet().index(level0Element);
+ }
+
return img[eIdx] == marked;
};
// create the grid
this->gridPtr() = this->createGrid_(this->hostGrid_(), elementSelector, paramGroup);
+
this->loadBalance();
}
};
diff --git a/dumux/io/grid/gridmanager_yasp.hh b/dumux/io/grid/gridmanager_yasp.hh
index 6ecb00067e5b62781070e78a5bb0f70abf07374b..d98c73d279aebb191e04726c5fd07fe66c1915c8 100644
--- a/dumux/io/grid/gridmanager_yasp.hh
+++ b/dumux/io/grid/gridmanager_yasp.hh
@@ -94,10 +94,13 @@ public:
// get the overlap
const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1);
- // make the grid
- ParentType::gridPtr() = createGrid_(modelParamGroup, upperRight, cells, periodic, overlap, Coordinates{});
-
- postProcessing_(modelParamGroup);
+ if constexpr (std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>)
+ init(upperRight, cells, modelParamGroup, overlap, periodic);
+ else
+ {
+ const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0));
+ init(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic);
+ }
}
// Didn't find a way to construct the grid
@@ -111,57 +114,67 @@ public:
}
}
-private:
/*!
- * \brief Create a grid with zero offset
+ * \brief Make the grid using input data not read from the input file.
+ * \note Use this function for EquidistantCoordinates where lowerLeft is always zero.
*/
- std::unique_ptr<Grid> createGrid_(const std::string& modelParamGroup,
- const GlobalPosition& upperRight,
- const std::array<int, dim>& cells,
- const std::bitset<dim>& periodic,
- const int overlap,
- Dune::EquidistantCoordinates<ct, dim>) const
+ void init(const GlobalPosition& upperRight,
+ const std::array<int, dim>& cells,
+ const std::string& modelParamGroup = "",
+ const int overlap = 1,
+ const std::bitset<dim> periodic = std::bitset<dim>{})
{
+ static_assert(std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>,
+ "Use init function taking lowerLeft as argument when working with EquidistantOffsetCoordinates");
+
if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning"))
{
// construct using default load balancing
- return std::make_unique<Grid>(upperRight, cells, periodic, overlap);
+ ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap);
}
else
{
// construct using user defined partitioning
const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
Dune::YaspFixedSizePartitioner<dim> lb(partitioning);
- return std::make_unique<Grid>(upperRight, cells, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
+ ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
}
+
+ postProcessing_(modelParamGroup);
}
/*!
- * \brief Create a grid with non-zero offset
+ * \brief Make the grid using input data not read from the input file.
+ * \note Use this function for EquidistantOffsetCoordinates.
*/
- std::unique_ptr<Grid> createGrid_(const std::string& modelParamGroup,
- const GlobalPosition& upperRight,
- const std::array<int, dim>& cells,
- const std::bitset<dim>& periodic,
- const int overlap,
- Dune::EquidistantOffsetCoordinates<ct, dim>) const
+ void init(const GlobalPosition& lowerLeft,
+ const GlobalPosition& upperRight,
+ const std::array<int, dim>& cells,
+ const std::string& modelParamGroup = "",
+ const int overlap = 1,
+ const std::bitset<dim> periodic = std::bitset<dim>{})
{
- const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0));
+ static_assert(std::is_same_v<Dune::EquidistantOffsetCoordinates<ct, dim>, Coordinates>,
+ "LowerLeft can only be specified with EquidistantOffsetCoordinates");
if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning"))
{
// construct using default load balancing
- return std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap);
+ ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap);
}
else
{
// construct using user defined partitioning
const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
Dune::YaspFixedSizePartitioner<dim> lb(partitioning);
- return std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
+ ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
}
+
+ postProcessing_(modelParamGroup);
}
+private:
+
/*!
* \brief Postprocessing for YaspGrid
*/
@@ -252,9 +265,7 @@ public:
const std::array<std::vector<ctype>, dim>& grading,
const std::string& modelParamGroup = "")
{
-
-
- // Additional arameters (they have a default)
+ // Additional parameters (they have a default)
const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1);
const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
// \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid)