diff --git a/src/sand-wedge-data/mygrid.cc b/src/sand-wedge-data/mygrid.cc
index 8f41ddb8343e214de92e27722f7a31148c922c33..d4fa12eac3ab5968ef80af8c7a49c55400d4e2a7 100644
--- a/src/sand-wedge-data/mygrid.cc
+++ b/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>
diff --git a/src/sand-wedge-data/mygrid.hh b/src/sand-wedge-data/mygrid.hh
index 1f0c5688677da3195e73b5d56249fd07d9680cee..cd17d0858f0475ac736f578eabc663a3bd79f20f 100644
--- a/src/sand-wedge-data/mygrid.hh
+++ b/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;
diff --git a/src/sand-wedge-data/mygrid_tmpl.cc b/src/sand-wedge-data/mygrid_tmpl.cc
index 3afed427fc51b0dd542d6a01e316fc4b5ce35d72..2668596c42fe202deeb150408066bc48508f7e66 100644
--- a/src/sand-wedge-data/mygrid_tmpl.cc
+++ b/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>;
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 6efacbc27a0bf00592d97a64e42fecc687aa2004..50326584720c881b0116bc713e44bd0eaf4d62e1 100644
--- a/src/sand-wedge.cc
+++ b/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);