[test] Add test for creating subgrid from binary image

Based on initial work by Kilian Weishaupt

test/io/gridmanager/CMakeLists.txt

@@ -35,7 +35,9 @@ dune_add_test(NAME test_gridmanager_subgrid
               CMD_ARGS --script fuzzy
                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_gridmanager_subgrid"
                        --files ${CMAKE_SOURCE_DIR}/test/references/subgrid-reference.vtu
-                               ${CMAKE_CURRENT_BINARY_DIR}/subgrid_three.vtu)
+                               ${CMAKE_CURRENT_BINARY_DIR}/subgrid_three.vtu
+                               ${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_subgrid_binary_image.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/subgrid_binary_image.vtu)
 
 dune_symlink_to_source_files(FILES "test_gridmanager_gmsh_3d.input"
                                    "test_gridmanager_gmsh_e_markers.input"

test/io/gridmanager/grids/dumux_binary.pbm

+P1
+# CREATOR: GIMP PNM Filter Version 1.1
+100 25
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000111000100000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000001011001000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000011100000000000000000000000000000000000000000000000100000
+0000000000000000000000000000111000000000000010000000000000000000000000
+0000000000000000000000011000000000000000000000000100000001110000000000
+0001100000000000000000000000000000000000000000000011111111100000000000
+0000011111000000110000000000000010110000000000000000000000000000000000
+0000000001101110011100000000000000111111000001110000000000000100110000
+0000000000000000000000000000000000000011001100001110000000000001100111
+0000111100000000000011001100000000000000000000000000000000000000000110
+0011000011100011001100110001110000111000001100110000000100000000000000
+0000000000000000000000000001100011000011100111001100110001111001011000
+0111001100000001100000000000000000000000000000000000000001100111000011
+1011100110001110011110111110001110011000000000111100000000000000000000
+0000000000000000010001100000110001100110000100110111101100000110011000
+0000000000000000000000000000000000000000000000000001100000110001100110
+0000001101110011000001100110000000000000000000000000000000000000000000
+0000000000011000011000011011100000001001110011000001101110000000000000
+0000000000000000000000000000000000000000110000111000111111111000011001
+1001110000111111111000000000000000000000000000000000000000000000000011
+1111111000001111111100111100011001111110111111110000000000000000000000
+0000000000000000000000000001111111100000001100110001111000010001111000
+1100110000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000000000
+00000000000000000000000000000000000000000000000000

test/io/gridmanager/test_gridmanager_subgrid.cc

@@ -20,6 +20,7 @@
  */
 #include <config.h>
 #include <iostream>
+#include <cmath>
 
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/fvector.hh>
@@ -48,7 +49,7 @@ public:
         const auto x = element.geometry().center()[0];
         const auto y = element.geometry().center()[1];
         const double radius = 0.3;
-        return std::sqrt((x - center_[0])*(x - center_[0]) + (y - center_[1])*(y - center_[1])) < radius;
+        return std::hypot(x-center_[0], y-center_[1]) < radius;
     }
 private:
     const GlobalPosition center_;
@@ -71,9 +72,9 @@ int main(int argc, char** argv) try
     using HostGrid = Dune::YaspGrid<dim, Dune::TensorProductCoordinates<double, dim> >;
 
     using HostGridManager = Dumux::GridManager<HostGrid>;
-    HostGridManager hostGridManager;
-    hostGridManager.init();
-    auto& hostGrid = hostGridManager.grid();
+    HostGridManager externalHostGridManager;
+    externalHostGridManager.init("External");
+    auto& hostGrid = externalHostGridManager.grid();
 
     // Calculate the bounding box of the host grid view.
     GlobalPosition bBoxMin(std::numeric_limits<double>::max());
@@ -92,58 +93,78 @@ int main(int argc, char** argv) try
     // Get the center of the hostgrid's domain.
     const GlobalPosition center{bBoxMin[0]+0.5*bBoxMax[0], bBoxMin[1]+0.5*bBoxMax[1]};
 
-    // Select all elements right of the center.
-    auto elementSelectorOne = [&center](const auto& element)
-    {
-        return element.geometry().center()[0] > center[0];
-    };
-
-    // Select all elements left of the center.
-    auto elementSelectorTwo = [&center](const auto& element)
-    {
-        return element.geometry().center()[0] < center[0];
-    };
-
-    // Select all elements within a circle around the center.
-    // Instead of a lambda, we use a class providing an () operator.
-    // Of course, a lambda would be possible here, too.
-    CircleSelector<GlobalPosition> elementSelectorThree(center);
-
-    // Create three different subgrids from the same hostgrid.
-    auto subgridPtrOne = SubgridManager<HostGrid>::makeGrid(hostGrid, elementSelectorOne);
-    auto subgridPtrTwo = SubgridManager<HostGrid>::makeGrid(hostGrid, elementSelectorTwo);
-    auto subgridPtrThree = SubgridManager<HostGrid>::makeGrid(hostGrid, elementSelectorThree);
-
-    std::cout << "Constructing a host grid and three subgrids took "  << timer.elapsed() << " seconds.\n";
-
     // Write out the host grid and the subgrids.
     {
         Dune::VTKWriter<HostGrid::LeafGridView> vtkWriter(hostGrid.leafGridView());
         vtkWriter.write("hostgrid");
     }
+
+    // Create different subgrids from the same hostgrid
     {
-        Dune::VTKWriter<SubgridManager<HostGrid>::Grid::LeafGridView> vtkWriter(subgridPtrOne->leafGridView());
+        std::cout << "Constructing SubGrid with tensor host grid and lambda element selector" << std::endl;
+
+        // Select all elements right of the center.
+        auto elementSelector = [&center](const auto& element)
+        {
+            return element.geometry().center()[0] > center[0];
+        };
+
+        SubgridManager<HostGrid> subgridManager;
+        subgridManager.init(hostGrid, elementSelector);
+        Dune::VTKWriter<SubgridManager<HostGrid>::Grid::LeafGridView> vtkWriter(subgridManager.grid().leafGridView());
         vtkWriter.write("subgrid_one");
     }
     {
-        Dune::VTKWriter<SubgridManager<HostGrid>::Grid::LeafGridView> vtkWriter(subgridPtrTwo->leafGridView());
+        std::cout << "Constructing SubGrid with tensor host grid and functor selector" << std::endl;
+
+        // Select all elements within a circle around the center.
+        // Instead of a lambda, we use a class providing an () operator.
+        // Of course, a lambda would be possible here, too.
+        CircleSelector<GlobalPosition> elementSelector(center);
+
+        SubgridManager<HostGrid> subgridManager;
+        subgridManager.init(hostGrid, elementSelector);
+        Dune::VTKWriter<SubgridManager<HostGrid>::Grid::LeafGridView> vtkWriter(subgridManager.grid().leafGridView());
+        vtkWriter.write("subgrid_three");
+    }
+
+    // create without contructing host grid first
+    {
+        std::cout << "Constructing SubGrid from lambda without specifying host grid" << std::endl;
+
+        // Select all elements left of the center.
+        auto elementSelector = [&center](const auto& element)
+        {
+            return element.geometry().center()[0] < center[0];
+        };
+
+        SubgridManager<HostGrid, HostGridManager> subgridManager;
+        subgridManager.init(elementSelector, "Internal");
+        Dune::VTKWriter<SubgridManager<HostGrid>::Grid::LeafGridView> vtkWriter(subgridManager.grid().leafGridView());
         vtkWriter.write("subgrid_two");
     }
+
+    // create subgrid from image file
     {
-        Dune::VTKWriter<SubgridManager<HostGrid>::Grid::LeafGridView> vtkWriter(subgridPtrThree->leafGridView());
-        vtkWriter.write("subgrid_three");
+        std::cout << "Constructing SubGrid from binary image" << std::endl;
+        using GridManager = SubgridManager<Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>>;
+        GridManager subgridManager; subgridManager.init("Image");
+        Dune::VTKWriter<GridManager::Grid::LeafGridView> vtkWriter(subgridManager.grid().leafGridView());
+        vtkWriter.write("subgrid_binary_image");
     }
 
+    std::cout << "Constructing a host grid and four subgrids took "  << timer.elapsed() << " seconds.\n";
+
     return 0;
 }
 ///////////////////////////////////////
 //////// Error handler ////////////////
 ///////////////////////////////////////
-catch (Dumux::ParameterException &e) {
+catch (const Dumux::ParameterException& e) {
     std::cerr << e << ". Abort!\n";
     return 1;
 }
-catch (Dune::Exception &e) {
+catch (const Dune::Exception& e) {
     std::cerr << "Dune reported error: " << e << std::endl;
     return 3;
 }

test/io/gridmanager/test_gridmanager_subgrid.input

-SubGridOne.Problem.Name = one
-SubGridTwo.Problem.Name = two
-SubGridThree.Problem.Name = three
+[External.Grid]
+Positions0 = 0 1
+Positions1 = 0 1
+Cells0 = 1
+Cells1 = 1
+Refinement = 4
 
-[Grid]
+[Internal.Grid]
 Positions0 = 0 1
 Positions1 = 0 1
 Cells0 = 1
 Cells1 = 1
 Refinement = 4
+
+[Image.Grid]
+LowerLeft = 1 1
+UpperRight = 5 2
+Cells = 100 25
+Image = grids/dumux_binary.pbm