From e593251eda96303b6f20d8fbc9a121b46a6b8202 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 25 Jun 2014 21:00:55 +0200
Subject: [PATCH] [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 | 78 ++++++++++++++++++-----------------
src/sand-wedge-data/mygrid.hh | 10 ++++-
2 files changed, 49 insertions(+), 39 deletions(-)
diff --git a/src/sand-wedge-data/mygrid.cc b/src/sand-wedge-data/mygrid.cc
index 2ef4059b..8f41ddb8 100644
--- a/src/sand-wedge-data/mygrid.cc
+++ b/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"
diff --git a/src/sand-wedge-data/mygrid.hh b/src/sand-wedge-data/mygrid.hh
index 9cdde76a..1f0c5688 100644
--- a/src/sand-wedge-data/mygrid.hh
+++ b/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
--
GitLab