[Extend] Add 3D sand wedge problem

.gitignore

@@ -26,5 +26,5 @@
 /test-driver
 Makefile
 Makefile.in
-src/sand-wedge
+src/sand-wedge-?D
 src/sliding-block-?D

src/Makefile.am

-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 = \

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);

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) {

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:

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