Grid refinement / adaptivity broken for level > 1
I've used subgrid to cut out some elements from a hostgrid, based on the hostgrid element indices (2D Yasp hostgrid).
Afterwards, I wanted to globally refine the subgrid.
I've tried
const auto numRefinements = 1 // 2, 3 ....
for (int i = 0; i < numRefinements; ++i)
{
for (const auto& e : elements(grid.leafGridView()))
grid.mark(1, e);
grid.preAdapt();
grid.adapt();
grid.postAdapt();
}
which is similiar to what is done in one of the tests.
As an alternative, I've also used
const auto numRefinements = 1 // 2, 3 ....
grid.globalRefine(numRefinements);
Both options work fine for numRefinements == 1
but fail for multiple refinements (numRefinements > 1
).
The second approach always seems to refine only once, while the first one does not preserve the original geometry of the cut-out regions.
Am I mis-using the refiment mechanism or could this be a bug?
Edit
Here is a MWE:
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <iostream>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/common/exceptions.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/subgrid/subgrid.hh>
using namespace Dune;
auto testSubGrid(const std::string& refinementStrategy)
{
using HostGrid = Dune::YaspGrid<2>;
using SubGrid = Dune::SubGrid<2, HostGrid>;
auto hostGrid = StructuredGridFactory<HostGrid>::createCubeGrid({0,0},{1,1},{3,3});
auto subGrid = SubGrid(*hostGrid);
//
// // A container to store the host grid elements' ids.
std::set<typename HostGrid::Traits::GlobalIdSet::IdType> elementsForSubgrid;
const auto& globalIDset = subGrid.getHostGrid().globalIdSet();
// Construct the subgrid.
subGrid.createBegin();
auto hostGridView = subGrid.getHostGrid().leafGridView();
const auto selector = [&](const auto& element)
{
// cut out element in center of 3x3 domain
return hostGridView.indexSet().index(element) != 4;
};
// Loop over all elements of the host grid and use the selector to
// choose which elements to add to the subgrid.
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");
subGrid.insertSetPartial(elementsForSubgrid);
subGrid.createEnd();
Dune::VTKWriter<typename SubGrid::LeafGridView> vtkWriter(subGrid.leafGridView());
vtkWriter.write("subgrid_unrefined");
if (refinementStrategy == "globalRefine")
{
subGrid.globalRefine(1);
vtkWriter.write("subgrid_1x_refined_globalRefine");
subGrid.globalRefine(1);
vtkWriter.write("subgrid_2x_refined_globalRefine");
}
if (refinementStrategy == "loop")
{
for (int i = 0; i < 2; ++i)
{
for (const auto& e : elements(subGrid.leafGridView()))
subGrid.mark(1, e);
subGrid.preAdapt();
subGrid.adapt();
subGrid.postAdapt();
vtkWriter.write("subgrid_" + std::to_string(i+1) + "x_refined_loop");
}
}
}
int main (int argc , char **argv) try
{
Dune::MPIHelper::instance(argc, argv);
testSubGrid("globalRefine");
testSubGrid("loop");
return 0;
}
catch (Dune::Exception& e) {
std::cerr << e << std::endl;
return 1;
} catch (...) {
std::cerr << "Generic exception!" << std::endl;
return 2;
}
Edited by Timo Koch