[Cleanup] Introduce SimplexManager class

src/sand-wedge-data/mygrid.cc

@@ -3,60 +3,122 @@
 #endif
 
 #include "mygrid.hh"
+#include "midpoint.hh"
 
-template <class Grid> std::shared_ptr<Grid> constructGrid() {
-  Dune::GridFactory<Grid> gridFactory;
+class SimplexManager {
+public:
+  using SimplexList = std::vector<std::vector<unsigned int>>;
 
 #if DIM == 3
-  Dune::FieldMatrix<double, 6, DIM> vertices;
-#else
-  Dune::FieldMatrix<double, 3, DIM> vertices;
+  SimplexManager(unsigned int shift) : shift_(shift) {}
 #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];
+
+  void addFromVertices(unsigned int U, unsigned int V, unsigned int W) {
 #if DIM == 3
-    vertices[3][i] = MyGeometry::A[i];
-    vertices[4][i] = MyGeometry::B[i];
-    vertices[5][i] = MyGeometry::C[i];
+    unsigned int const U2 = U + shift_;
+    unsigned int const V2 = V + shift_;
+    unsigned int const W2 = W + shift_;
+
+    // back-to-front, back-to-front, back-to-front
+    simplices_.push_back({ U, V, W, U2 });
+    simplices_.push_back({ V, V2, W2, U2 });
+    simplices_.push_back({ W, W2, U2, V });
+
+// back-to-front, front-to-back, back-to-front would be
+/*
+  simplices_.push_back({ U, V, W, U2 });
+  simplices_.push_back({ V, V2, W, U2 });
+  simplices_.push_back({ V2, W2, U2, W });
+*/
+#else
+    simplices_.push_back({ U, V, W });
 #endif
   }
 
+  SimplexList const &getSimplices() { return simplices_; }
+
+private:
+  SimplexList simplices_;
+
 #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;
-  }
+  unsigned int const shift_;
+#endif
+};
+
+template <class Grid> class GridConstructor {
+public:
+  GridConstructor() {
+    auto const &A = MyGeometry::A;
+    auto const &B = MyGeometry::B;
+    auto const &C = MyGeometry::C;
+
+    unsigned int const vc = 3;
+
+#if DIM == 3
+    Dune::FieldMatrix<double, 2 * vc, DIM> vertices;
+#else
+    Dune::FieldMatrix<double, vc, DIM> vertices;
+#endif
+    for (size_t i = 0; i < 2; ++i) {
+      size_t k = 0;
+      vertices[k++][i] = A[i];
+      vertices[k++][i] = B[i];
+      vertices[k++][i] = C[i];
+      assert(k == vc);
+#if DIM == 3
+      vertices[k++][i] = A[i];
+      vertices[k++][i] = B[i];
+      vertices[k++][i] = C[i];
+      assert(k == 2 * vc);
+#endif
+    }
+
+#if DIM == 3
+    for (size_t k = 0; k < vc; ++k) {
+      vertices[k][2] = -MyGeometry::depth / 2.0;
+      vertices[k + vc][2] = MyGeometry::depth / 2.0;
+    }
 #endif
 
-  for (size_t i = 0; i < vertices.N(); ++i)
-    gridFactory.insertVertex(vertices[i]);
+    for (size_t i = 0; i < vertices.N(); ++i)
+      gridFactory.insertVertex(vertices[i]);
 
-  Dune::GeometryType cell;
+    Dune::GeometryType cell;
 #if DIM == 3
-  cell.makeTetrahedron();
+    cell.makeTetrahedron();
 #else
-  cell.makeTriangle();
+    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 } };
+    SimplexManager sm(vc);
 #else
-  std::vector<std::vector<unsigned int>> simplices = { { 0, 1, 2 } };
+    SimplexManager sm;
 #endif
+    sm.addFromVertices(0, 1, 2);
+    auto const &simplices = sm.getSimplices();
+
+    // sanity-check choices of simplices
+    for (size_t i = 0; i < simplices.size(); ++i) {
+      Dune::FieldMatrix<double, DIM, DIM> check;
+      for (size_t j = 0; j < DIM; ++j)
+        check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+      assert(check.determinant() > 0);
+      gridFactory.insertElement(cell, simplices[i]);
+    }
+  }
+
+  std::shared_ptr<Grid> getGrid() {
+    return std::shared_ptr<Grid>(gridFactory.createGrid());
+  }
 
-  // sanity-check choices of simplices
-  for (size_t i = 0; i < simplices.size(); ++i) {
-    Dune::FieldMatrix<double, DIM, DIM> check;
-    for (size_t j = 0; j < DIM; ++j)
-      check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
-    assert(check.determinant() > 0);
-    gridFactory.insertElement(cell, simplices[i]);
+  template <class GridView>
+  MyFaces<GridView> constructFaces(GridView const &gridView) {
+    return MyFaces<GridView>(gridView);
   }
-  return std::shared_ptr<Grid>(gridFactory.createGrid());
+
+private:
+  Dune::GridFactory<Grid> gridFactory;
 };
 
 template <class GridView>

src/sand-wedge-data/mygrid.hh

@@ -14,8 +14,6 @@
 
 #include "mygeometry.hh"
 
-template <class Grid> std::shared_ptr<Grid> constructGrid();
-
 template <class GridView> struct MyFaces {
   BoundaryPatch<GridView> lower;
   BoundaryPatch<GridView> right;

src/sand-wedge-data/mygrid_tmpl.cc

@@ -4,6 +4,6 @@
 
 #include "explicitgrid.hh"
 
-template std::shared_ptr<Grid> constructGrid();
+template class GridConstructor<Grid>;
 
 template struct MyFaces<GridView>;

src/sand-wedge.cc

@@ -106,7 +106,8 @@ int main(int argc, char *argv[]) {
 
     // {{{ Set up grid
     using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
-    auto grid = constructGrid<Grid>();
+    GridConstructor<Grid> gridConstructor;
+    auto grid = gridConstructor.getGrid();
 
     auto const refinements = parset.get<size_t>("grid.refinements");
     grid->globalRefine(refinements);
@@ -115,7 +116,7 @@ int main(int argc, char *argv[]) {
     using GridView = Grid::LeafGridView;
     GridView const leafView = grid->leafView();
 
-    MyFaces<GridView> myFaces(leafView);
+    auto myFaces = gridConstructor.constructFaces(leafView);
 
     // Neumann boundary
     BoundaryPatch<GridView> const neumannBoundary(leafView);