From d0deb52798fbe23e77019eca54cf83d6fc7d24a3 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 17 Jun 2014 14:14:02 +0200
Subject: [PATCH] [Extend] Add 3D sand wedge problem
---
.gitignore | 2 +-
src/Makefile.am | 10 +++--
src/sand-wedge-data/mygeometry.hh | 2 +
src/sand-wedge-data/mygrid.cc | 61 +++++++++++++++++++++++++++++--
src/sand-wedge-data/mygrid.hh | 6 +++
src/sand-wedge.cc | 5 +++
6 files changed, 78 insertions(+), 8 deletions(-)
diff --git a/.gitignore b/.gitignore
index f47390d2..10191e28 100644
--- a/.gitignore
+++ b/.gitignore
@@ -26,5 +26,5 @@
/test-driver
Makefile
Makefile.in
-src/sand-wedge
+src/sand-wedge-?D
src/sliding-block-?D
diff --git a/src/Makefile.am b/src/Makefile.am
index a939e404..45d92485 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = sand-wedge
+bin_PROGRAMS = sand-wedge-2D sand-wedge-3D
common_sources = \
assemblers.cc \
@@ -8,10 +8,14 @@ common_sources = \
timestepping.cc \
vtk.cc
-sand_wedge_SOURCES = $(common_sources) sand-wedge.cc
-sand_wedge_CPPFLAGS = \
+sand_wedge_2D_SOURCES = $(common_sources) sand-wedge.cc
+sand_wedge_2D_CPPFLAGS = \
$(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
-Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DDIM=2
+sand_wedge_3D_SOURCES = $(common_sources) sand-wedge.cc
+sand_wedge_3D_CPPFLAGS = \
+ $(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
+ -Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DDIM=3
# Some are for clang, others are for gcc
AM_CXXFLAGS = \
diff --git a/src/sand-wedge-data/mygeometry.hh b/src/sand-wedge-data/mygeometry.hh
index bf3baac5..468175d7 100644
--- a/src/sand-wedge-data/mygeometry.hh
+++ b/src/sand-wedge-data/mygeometry.hh
@@ -92,6 +92,8 @@ namespace {
}
}
+double const depth = 0.10;
+
LocalVector const A = rotate(reference::A);
LocalVector const B = rotate(reference::B);
LocalVector const C = rotate(reference::C);
diff --git a/src/sand-wedge-data/mygrid.cc b/src/sand-wedge-data/mygrid.cc
index ec234395..2ef4059b 100644
--- a/src/sand-wedge-data/mygrid.cc
+++ b/src/sand-wedge-data/mygrid.cc
@@ -7,17 +7,46 @@
template <class Grid> std::shared_ptr<Grid> constructGrid() {
Dune::GridFactory<Grid> gridFactory;
+#if DIM == 3
+ Dune::FieldMatrix<double, 6, DIM> vertices;
+#else
Dune::FieldMatrix<double, 3, DIM> vertices;
- vertices[0] = MyGeometry::A;
- vertices[1] = MyGeometry::B;
- vertices[2] = MyGeometry::C;
+#endif
+ for (size_t i = 0; i < 2; ++i) {
+ vertices[0][i] = MyGeometry::A[i];
+ vertices[1][i] = MyGeometry::B[i];
+ vertices[2][i] = MyGeometry::C[i];
+#if DIM == 3
+ vertices[3][i] = MyGeometry::A[i];
+ vertices[4][i] = MyGeometry::B[i];
+ vertices[5][i] = MyGeometry::C[i];
+#endif
+ }
+
+#if DIM == 3
+ for (size_t k = 0; k < 3; ++k) {
+ vertices[k][2] = -MyGeometry::depth / 2.0;
+ vertices[k + 3][2] = MyGeometry::depth / 2.0;
+ }
+#endif
+
for (size_t i = 0; i < vertices.N(); ++i)
gridFactory.insertVertex(vertices[i]);
Dune::GeometryType cell;
+#if DIM == 3
+ cell.makeTetrahedron();
+#else
cell.makeTriangle();
+#endif
+#if DIM == 3
+ std::vector<std::vector<unsigned int>> simplices = { { 0, 1, 2, 3 },
+ { 1, 2, 3, 4 },
+ { 2, 3, 4, 5 } };
+#else
std::vector<std::vector<unsigned int>> simplices = { { 0, 1, 2 } };
+#endif
// sanity-check choices of simplices
for (size_t i = 0; i < simplices.size(); ++i) {
@@ -32,16 +61,40 @@ template <class Grid> std::shared_ptr<Grid> constructGrid() {
template <class GridView>
MyFaces<GridView>::MyFaces(GridView const &gridView)
- : lower(gridView), right(gridView) {
+ :
+#if DIM == 3
+ lower(gridView),
+ right(gridView),
+ front(gridView),
+ back(gridView)
+#else
+ lower(gridView),
+ right(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) {
return isClose(0, in.geometry().center()[1]);
});
+#if DIM == 3
+ front.insertFacesByProperty([&](
+ typename Grid::LeafGridView::Intersection const &in) {
+ return isClose(-MyGeometry::depth / 2.0, in.geometry().center()[2]);
+ });
+
+ back.insertFacesByProperty([&](
+ typename Grid::LeafGridView::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) {
diff --git a/src/sand-wedge-data/mygrid.hh b/src/sand-wedge-data/mygrid.hh
index 1f02943d..9cdde76a 100644
--- a/src/sand-wedge-data/mygrid.hh
+++ b/src/sand-wedge-data/mygrid.hh
@@ -21,6 +21,12 @@ template <class GridView> struct MyFaces {
BoundaryPatch<GridView> lower;
BoundaryPatch<GridView> right;
BoundaryPatch<GridView> upper;
+
+#if DIM == 3
+ BoundaryPatch<GridView> front;
+ BoundaryPatch<GridView> back;
+#endif
+
MyFaces(GridView const &gridView);
private:
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 264eb169..c58817d7 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -140,6 +140,11 @@ int main(int argc, char *argv[]) {
if (myFaces.lower.containsVertex(i))
dirichletNodes[i][1] = true;
+
+#if DIM == 3
+ if (myFaces.front.containsVertex(i) || myFaces.back.containsVertex(i))
+ dirichletNodes[i][2] = true;
+#endif
}
// Set up functions for time-dependent boundary conditions
--
GitLab