[Extend] Determine boundary sections properly

The previous solution assumed that some sections consisted of a single
face on the coarse grid

src/sand-wedge-data/mygrid.cc

@@ -57,6 +57,35 @@ template <class Grid> std::shared_ptr<Grid> constructGrid() {
     gridFactory.insertElement(cell, simplices[i]);
   }
   return std::shared_ptr<Grid>(gridFactory.createGrid());
+};
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+                                    Vector const &c) {
+  return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+                                Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 > x[0] or x[0] > minmax0.second + 1e-14)
+    return false;
+  if (minmax1.first - 1e-14 > x[1] or x[1] > minmax1.second + 1e-14)
+    return false;
+
+  return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+                                  Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
 }
 
 template <class GridView>
@@ -65,65 +94,38 @@ MyFaces<GridView>::MyFaces(GridView const &gridView)
 #if DIM == 3
       lower(gridView),
       right(gridView),
+      upper(gridView),
       front(gridView),
       back(gridView)
 #else
       lower(gridView),
-      right(gridView)
+      right(gridView),
+      upper(gridView)
 #endif
 {
   using Grid = typename GridView::Grid;
-  using LevelGridView = typename Grid::LevelGridView;
-  LevelGridView const rootView = gridView.grid().levelView(0);
-
   assert(isClose(MyGeometry::A[1], 0));
   assert(isClose(MyGeometry::B[1], 0));
-  lower.insertFacesByProperty([&](
-      typename Grid::LeafGridView::Intersection const &in) {
+  lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
     return isClose(0, in.geometry().center()[1]);
   });
-
 #if DIM == 3
-  front.insertFacesByProperty([&](
-      typename Grid::LeafGridView::Intersection const &in) {
+  front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
     return isClose(-MyGeometry::depth / 2.0, in.geometry().center()[2]);
   });
 
-  back.insertFacesByProperty([&](
-      typename Grid::LeafGridView::Intersection const &in) {
+  back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
     return isClose(MyGeometry::depth / 2.0, in.geometry().center()[2]);
   });
 #endif
 
-  BoundaryPatch<LevelGridView> upperRoot(rootView);
-  upperRoot.insertFacesByProperty([&](
-      typename LevelGridView::Intersection const &in) {
-    auto const geometry = in.geometry();
-    auto const range = boost::irange(0, geometry.corners());
-    return boost::algorithm::all_of(range, [&](int i) {
-      auto const &co = geometry.corner(i);
-      return (isClose(co[1], MyGeometry::A[1]) &&
-              isClose(co[0], MyGeometry::A[0])) ||
-             (isClose(co[1], MyGeometry::C[1]) &&
-              isClose(co[0], MyGeometry::C[0]));
-    });
+  upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return xyBetween(MyGeometry::A, MyGeometry::C, in.geometry().center());
   });
-  BoundaryPatchProlongator<Grid>::prolong(upperRoot, upper);
-
-  BoundaryPatch<LevelGridView> rightRoot(rootView);
-  rightRoot.insertFacesByProperty([&](
-      typename LevelGridView::Intersection const &in) {
-    auto const geometry = in.geometry();
-    auto const range = boost::irange(0, geometry.corners());
-    return boost::algorithm::all_of(range, [&](int i) {
-      auto const &co = geometry.corner(i);
-      return (isClose(co[1], MyGeometry::C[1]) &&
-              isClose(co[0], MyGeometry::C[0])) ||
-             (isClose(co[1], MyGeometry::B[1]) &&
-              isClose(co[0], MyGeometry::B[0]));
-    });
+
+  right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return xyBetween(MyGeometry::B, MyGeometry::C, in.geometry().center());
   });
-  BoundaryPatchProlongator<Grid>::prolong(rightRoot, right);
 }
 
 #include "mygrid_tmpl.cc"

src/sand-wedge-data/mygrid.hh

@@ -11,7 +11,6 @@
 #pragma clang diagnostic ignored "-Wdeprecated-declarations"
 #include <dune/fufem/boundarypatch.hh>
 #pragma clang diagnostic pop
-#include <dune/fufem/boundarypatchprolongator.hh>
 
 #include "mygeometry.hh"
 
@@ -33,5 +32,14 @@ template <class GridView> struct MyFaces {
   bool isClose(double a, double b) {
     return std::abs(a - b) < 1e-14;
   };
+
+  template <class Vector>
+  bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
+
+  template <class Vector>
+  bool xyCollinear(Vector const &a, Vector const &b, Vector const &c);
+
+  template <class Vector>
+  bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
 };
 #endif