Skip to content

Grid refinement seems to be 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;
}