diff --git a/dune/tectonic/problem-data/grid/mygrids.cc b/dune/tectonic/problem-data/grid/mygrids.cc index e3179f2a84c78e6eb87bd920c9c458ab68261956..4e0210f9bb4ea86c081f04f307586a3833229b21 100644 --- a/dune/tectonic/problem-data/grid/mygrids.cc +++ b/dune/tectonic/problem-data/grid/mygrids.cc @@ -10,9 +10,62 @@ #include "simplexmanager.hh" + + + + +// ------------------------------------------------- +// GridsConstructor +// ------------------------------------------------- + +template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries) : + cuboidGeometries_(cuboidGeometries) +{ + size_t const gridCount = cuboidGeometries_.size(); + grids_.resize(gridCount); + gridFactories_.resize(gridCount); + + buildGrids(); +} + +template <class Grid> +std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() { + return grids; +} + +template <class Grid> +template <class GridView> +MyFaces<GridView> GridsConstructor<Grid>::constructFaces( + GridView const &gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry) { + return MyFaces<GridView>(gridView, cuboidGeometry); +} + + + // Fix: 3D case (still Elias code) template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_) : cuboidGeometries(cuboidGeometries_) + + +template <class Grid> +std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() { + return grids; +} + +template <class Grid> +template <class GridView> +MyFaces<GridView> GridsConstructor<Grid>::constructFaces( + GridView const &gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry) { + return MyFaces<GridView>(gridView, cuboidGeometry); +} + + +// ------------------------------------------------- +// UniformGridsContructor +// ------------------------------------------------- + +template <class Grid> UniformGridsConstructor<Grid>::UniformGridsContructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_) : + cuboidGeometries(cuboidGeometries_) { size_t const gridCount = cuboidGeometries.size(); grids.resize(gridCount); @@ -90,7 +143,7 @@ template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector } template <class Grid> -std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() { +std::vector<std::shared_ptr<Grid>>& UniformGridsContructor<Grid>::getGrids() { return grids; } @@ -101,81 +154,6 @@ MyFaces<GridView> GridsConstructor<Grid>::constructFaces( return MyFaces<GridView>(gridView, cuboidGeometry); } -template <class GridView> -template <class Vector> -bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b, - Vector const &c) { - return isClose2((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 * cuboidGeometry.lengthScale() > x[0] or - x[0] > minmax0.second + 1e-14 * cuboidGeometry.lengthScale()) - return false; - if (minmax1.first - 1e-14 * cuboidGeometry.lengthScale() > x[1] or - x[1] > minmax1.second + 1e-14 * cuboidGeometry.lengthScale()) - 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> -MyFaces<GridView>::MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_) - : - #if MY_DIM == 3 - lower(gridView), - right(gridView), - upper(gridView), - left(gridView), - front(gridView), - back(gridView), - #else - lower(gridView), - right(gridView), - upper(gridView), - left(gridView), - #endif - cuboidGeometry(cuboidGeometry_) - { - lower.insertFacesByProperty([&](typename GridView::Intersection const &in) { - return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.lowerRight(), in.geometry().center()); - }); - - right.insertFacesByProperty([&](typename GridView::Intersection const &in) { - return xyBetween(cuboidGeometry.lowerRight(), cuboidGeometry.upperRight(), in.geometry().center()); - }); - - upper.insertFacesByProperty([&](typename GridView::Intersection const &in) { - return xyBetween(cuboidGeometry.upperLeft(), cuboidGeometry.upperRight(), in.geometry().center()); - }); - - left.insertFacesByProperty([&](typename GridView::Intersection const &in) { - return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.upperLeft(), in.geometry().center()); - }); - - #if MY_DIM == 3 - front.insertFacesByProperty([&](typename GridView::Intersection const &in) { - return isClose(cuboidGeometry.depth() / 2.0, in.geometry().center()[2]); - }); - - back.insertFacesByProperty([&](typename GridView::Intersection const &in) { - return isClose(-cuboidGeometry.depth() / 2.0, in.geometry().center()[2]); - }); - #endif - } double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) { return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter; diff --git a/dune/tectonic/problem-data/grid/mygrids.hh b/dune/tectonic/problem-data/grid/mygrids.hh index 62b440f5068d51827bcd28065258df13b989550a..6883c228041c73e1a9b02c0013d23eb533e2ccbb 100644 --- a/dune/tectonic/problem-data/grid/mygrids.hh +++ b/dune/tectonic/problem-data/grid/mygrids.hh @@ -44,11 +44,11 @@ template <class GridView> struct MyFaces { bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x); }; -template <class Grid> class GridsConstructor { -public: - GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_); - std::vector<std::shared_ptr<Grid>>& getGrids(); +template <class Grid> +class CuboidGridConstructor : public GridConstructor<Grid> { +public: + CuboidGridConstructor(const CuboidGeometry<typename Grid::ctype>& cuboidGeometry); template <class GridView> MyFaces<GridView> constructFaces(const GridView& gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry); @@ -59,12 +59,23 @@ template <class Grid> class GridsConstructor { std::vector<std::shared_ptr<Grid>> grids; }; + +template <class Grid> +class UniformGridsConstructor : public GridsConstructor<Grid> { +public: + UniformGridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_); +}; + + double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale); template <class Grid, class LocalVector> void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch, double smallestDiameter, double lengthScale); +template <class Grid> +void refine(Grid &grid, size_t uniformRefinements); + #endif