Merge branch 'bugfix/issue-53-uggrid-intersections' into 'master'

UGGrid: append to-be-visited sons at the end of the list

Closes #53

See merge request !150

dune/grid/test/CMakeLists.txt

@@ -41,6 +41,9 @@ dune_add_test(SOURCES test-parallel-ug.cc
 dune_add_test(SOURCES test-loadbalancing.cc
               CMAKE_GUARD UG_FOUND)
 
+dune_add_test(SOURCES issue-53-uggrid-intersections.cc
+              CMAKE_GUARD UG_FOUND)
+
 # The alberta tests are only alibi-ported, until the  grid and world dimension
 # are configuretime parameters and we can treat alberta just as any other grid manager
 # - buildsystemwise. PLEASE DON'T LOOK AT THIS IF YOU WANT TO KNOW HOW TO WRITE TESTS.

dune/grid/test/issue-53-uggrid-intersections.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+/*
+ * Check leaf intersection geometry on UGGrid.
+ *
+ * Reference: https://gitlab.dune-project.org/core/dune-grid/issues/53
+ */
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+#include <iostream>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/vtk.hh>
+#if HAVE_DUNE_ALUGRID
+#  include <dune/alugrid/grid.hh>
+#endif
+
+const double EPSILON = 1e-12;
+
+using Dune::TestSuite;
+
+template<class Grid>
+TestSuite testGrid(const Grid& grid, const std::string& name) {
+  TestSuite test;
+
+  // now iterate over intersections
+  const auto& gv = leafGridView(grid);
+  for (const auto& e : elements(gv)) {
+    for (const auto& is : intersections(gv, e)) {
+      if (is.boundary())
+        continue;
+
+      // compute the edge lengths of the inner and outer element
+      auto insideLength = std::sqrt(is.inside().geometry().volume());
+      auto outsideLength = std::sqrt(is.outside().geometry().volume());
+      // length of intersection
+      auto intersectionLength = is.geometry().volume();
+
+      // intersection should not be longer than any full edge
+      const bool check = intersectionLength < insideLength+EPSILON && intersectionLength < outsideLength+EPSILON;
+      test.check(check)
+        << "=== " << name << ": Intersection is longer than actual edge! ===\n"
+        << "  inside center: " << is.inside().geometry().center() << "\n"
+        << " outside center: " << is.outside().geometry().center() << "\n"
+        << " lengths: intersection: " << intersectionLength << " inside: " << insideLength << " outside: " << outsideLength << "\n";
+    }
+  }
+
+  // output mesh as VTK (uncomment if needed)
+  // auto vtk = Dune::VTKWriter<typename Grid::LeafGridView>(gv);
+  // vtk.write(s);
+
+  return test;
+}
+
+template<class Grid>
+TestSuite testRefinement1(Grid& grid, const std::string& name)
+{
+  const auto& gv = leafGridView(grid);
+
+  // Make some pseudo-adaptive refinements
+  for (int i = 0; i < 5; ++i)
+  {
+    for (auto e : elements(gv)) {
+      if (gv.indexSet().subIndex(e,0,0) % 5 == 0) // mark some elements for refinement
+        grid.mark(1, e);
+    }
+    grid.adapt();
+    grid.postAdapt();
+  }
+
+  return testGrid(grid, name);
+}
+
+template<typename Grid>
+TestSuite testRefinement2(Grid& grid, const std::string& name)
+{
+  using Coordinate = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
+
+  // List of list of elements to refine in each refinement cycle.
+  // The elements are specified using their cell center.
+  const std::vector< std::vector<Coordinate> > refine = {
+    {{0.5, 0.5}},
+    {{0.75, 0.25}},
+    {{0.875, 0.375}},
+    {{0.9375, 0.4375}, {0.875, 0.125}}
+  };
+
+  const auto& gv = leafGridView(grid);
+
+  for (const auto& xs : refine) {
+    for (const auto& e : elements(gv)) {
+      for (const auto& x : xs) {
+        auto d = e.geometry().center() - x;
+        if (d.two_norm() < EPSILON)
+          grid.mark(1, e);
+      }
+    }
+    grid.adapt();
+    grid.postAdapt();
+  }
+
+  return testGrid(grid, name);
+}
+
+int main(int argc, char** argv)
+{
+  Dune::MPIHelper::instance(argc, argv);
+
+  TestSuite test;
+
+  // test with UG
+  {
+    using GridType = Dune::UGGrid<2>;
+    auto grid = Dune::StructuredGridFactory<GridType>::createCubeGrid({0,0}, {1.0,1.0}, {16,16});
+    grid->setClosureType(GridType::ClosureType::NONE);
+    test.subTest(testRefinement1(*grid, "UGGrid-1"));
+  }
+  {
+    using GridType = Dune::UGGrid<2>;
+    auto grid = Dune::StructuredGridFactory<GridType>::createCubeGrid({0,0}, {2.0,1.0}, {2,1});
+    grid->setClosureType(GridType::ClosureType::NONE);
+    test.subTest(testRefinement2(*grid, "UGGrid-2"));
+  }
+  // test with ALU
+#if HAVE_DUNE_ALUGRID
+  {
+    using GridType = Dune::ALUGrid<2,2, Dune::cube, Dune::nonconforming>;
+    auto grid = Dune::StructuredGridFactory<GridType>::createCubeGrid({0.0,0.0}, {1.0,1.0}, {16,16});
+    test.subTest(testRefinement1(*grid, "ALUGrid-1"));
+  }
+  {
+    using GridType = Dune::ALUGrid<2,2, Dune::cube, Dune::nonconforming>;
+    auto grid = Dune::StructuredGridFactory<GridType>::createCubeGrid({0,0}, {2.0,1.0}, {2,1});
+    grid->setClosureType(GridType::ClosureType::NONE);
+    test.subTest(testRefinement2(*grid, "ALUGrid-2"));
+  }
+#endif
+
+  return test.exit();
+}

dune/grid/uggrid/uggrid.cc

@@ -10,7 +10,7 @@
 #include <dune/grid/uggrid.hh>
 
 /** \todo Remove the following two includes once getAllSubfaces... is gone */
-#include <forward_list>
+#include <list>
 #include <iterator>
 #include <dune/common/stdstreams.hh>
 #include <dune/grid/common/mcmgmapper.hh>
@@ -360,7 +360,7 @@ void Dune::UGGrid<dim>::getChildrenOfSubface(const typename Traits::template Cod
 
   typedef std::pair<typename UG_NS<dim>::Element*,int> ListEntryType;
 
-  std::forward_list<ListEntryType> list;
+  std::list<ListEntryType> list;
 
   // //////////////////////////////////////////////////////////////////////
   //   Change the input face number from Dune numbering to UG numbering
@@ -375,23 +375,7 @@ void Dune::UGGrid<dim>::getChildrenOfSubface(const typename Traits::template Cod
       && e.level() < maxl) {
 
     typename UG_NS<dim>::Element* theElement = this->getRealImplementation(e).target_;
-
-    UG::INT Sons_of_Side = 0;
-    typename UG_NS<dim>::Element* SonList[UG_NS<dim>::MAX_SONS];
-    UG::INT SonSides[UG_NS<dim>::MAX_SONS];
-
-    int rv = Get_Sons_of_ElementSide(theElement,
-                                     elementSide,
-                                     &Sons_of_Side,
-                                     SonList,          // the output elements
-                                     SonSides,         // Output element side numbers
-                                     true,            // Element sons are not precomputed
-                                     true);            // ioflag: I have no idea what this is supposed to do
-    if (rv!=0)
-      DUNE_THROW(GridError, "Get_Sons_of_ElementSide returned with error value " << rv);
-
-    for (int i=0; i<Sons_of_Side; i++)
-      list.push_front(ListEntryType(SonList[i],SonSides[i]));
+    list.emplace_back(theElement, elementSide);
 
   }
 
@@ -410,21 +394,26 @@ void Dune::UGGrid<dim>::getChildrenOfSubface(const typename Traits::template Cod
 
     if (UG_NS<dim>::myLevel(theElement) < maxl) {
 
-      Get_Sons_of_ElementSide(theElement,
-                              side,         // Input element side number
-                              &Sons_of_Side,       // Number of topological sons of the element side
-                              SonList,            // Output elements
-                              SonSides,           // Output element side numbers
-                              true,
-                              true);
+      int rv = Get_Sons_of_ElementSide(theElement,
+                                       side,          // Input element side number
+                                       &Sons_of_Side, // Number of topological sons of the element side
+                                       SonList,       // Output elements
+                                       SonSides,      // Output element side numbers
+                                       true,
+                                       true);
+
+      if (rv != 0)
+        DUNE_THROW(GridError, "Get_Sons_of_ElementSide returned with error value " << rv);
 
       for (int i=0; i<Sons_of_Side; i++)
-        list.push_front(ListEntryType(SonList[i],SonSides[i]));
+        list.emplace_back(SonList[i], SonSides[i]);
 
     }
-
   }
 
+  // Remove the element itself. We are only interested in the children.
+  list.pop_front();
+
   // //////////////////////////////
   //   Extract result from stack
   // //////////////////////////////

dune/grid/uggrid/uggridintersections.cc

@@ -5,8 +5,8 @@
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/uggrid/uggridintersections.hh>
 
+#include <list>
 #include <set>
-#include <forward_list>
 
 template<class GridImp>
 const typename Dune::UGGridLevelIntersection<GridImp>::WorldVector&
@@ -577,7 +577,6 @@ int Dune::UGGridLeafIntersection<GridImp>::getFatherSide(const Face& currentFace
     }
 
     DUNE_THROW(InvalidStateException,"getFatherSide() didn't find a father.");
-    return 0;
 
   } else {    //  dim==3
 
@@ -634,7 +633,7 @@ int Dune::UGGridLeafIntersection<GridImp>::getFatherSide(const Face& currentFace
   }
 
   // Should never happen
-  return -1;
+  DUNE_THROW(GridError, "Reached code path that should never be reached");
 }
 
 template< class GridImp>
@@ -697,27 +696,10 @@ void Dune::UGGridLeafIntersection<GridImp>::constructLeafSubfaces() {
 
   else {
 
-    std::forward_list<Face> list;
-    int levelNeighborSide = numberInNeighbor(center_, levelNeighbor);
+    std::list<Face> list;
 
-    UG::INT Sons_of_Side = 0;
-    typename UG_NS<dim>::Element* SonList[UG_NS<dim>::MAX_SONS];
-    UG::INT SonSides[UG_NS<dim>::MAX_SONS];
-
-    int rv = Get_Sons_of_ElementSide(levelNeighbor,
-                                     levelNeighborSide,
-                                     &Sons_of_Side,
-                                     SonList,      // the output elements
-                                     SonSides,     // Output element side numbers
-                                     true,        // Element sons are not precomputed
-                                     false,        // ioflag: Obsolete debugging flag
-                                     true);
-
-    if (rv!=0)
-      DUNE_THROW(GridError, "Get_Sons_of_ElementSide returned with error value " << rv);
-
-    for (int i=0; i<Sons_of_Side; i++)
-      list.push_front(Face(SonList[i],SonSides[i]));
+    const int levelNeighborSide = numberInNeighbor(center_, levelNeighbor);
+    list.emplace_back(levelNeighbor, levelNeighborSide);
 
     // //////////////////////////////////////////////////
     //   Get_Sons_of_ElementSide only computes direct sons.  Therefore in order to get all
@@ -735,22 +717,30 @@ void Dune::UGGridLeafIntersection<GridImp>::constructLeafSubfaces() {
 
       if (!UG_NS<dim>::isLeaf(theElement)) {
 
-        Get_Sons_of_ElementSide(theElement,
-                                f.second,      // Input element side number
-                                &Sons_of_Side,     // Number of topological sons of the element side
-                                SonList,          // Output elements
-                                SonSides,         // Output element side numbers
-                                true,
-                                false,
-                                true);
+        int rv = Get_Sons_of_ElementSide(theElement,
+                                         f.second,      // Input element side number
+                                         &Sons_of_Side, // Number of topological sons of the element side
+                                         SonList,       // Output elements
+                                         SonSides,      // Output element side numbers
+                                         true,
+                                         false,
+                                         true);
+
+        if (rv!=0)
+          DUNE_THROW(GridError, "Get_Sons_of_ElementSide returned with error value " << rv);
 
         for (int i=0; i<Sons_of_Side; i++)
-          list.push_front(Face(SonList[i],SonSides[i]));
+          list.emplace_back(SonList[i], SonSides[i]);
 
       }
 
     }
 
+    // Remove the element itself.
+    // We are only interested in leaf elements and in this code branch
+    // the element we start with is never a leaf.
+    list.pop_front();
+
     // //////////////////////////////
     //   Extract result from stack
     // //////////////////////////////