[Cleanup] Print grid info

src/sand-wedge.cc

@@ -44,6 +44,7 @@
 #include <dune/fufem/sharedpointermap.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/norms/sumnorm.hh>
+#include <dune/solvers/norms/twonorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/solvers/solvers/solver.hh>
 #include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
@@ -81,6 +82,21 @@ void initPython() {
   Python::run("sys.path.append('" datadir "')");
 }
 
+template <class Geometry> double diameter(Geometry const &geometry) {
+  auto const numCorners = geometry.corners();
+  std::vector<typename Geometry::GlobalCoordinate> corners(numCorners);
+  for (int i = 0; i < numCorners; ++i)
+    corners[i] = geometry.corner(i);
+
+  TwoNorm<typename Geometry::GlobalCoordinate> twoNorm;
+
+  double diameter = 0.0;
+  for (int i = 0; i < numCorners; ++i)
+    for (int j = 0; j < i; ++j)
+      diameter = std::max(diameter, twoNorm.diff(corners[i], corners[j]));
+  return diameter;
+}
+
 int main(int argc, char *argv[]) {
   try {
     Dune::ParameterTree parset;
@@ -97,10 +113,24 @@ int main(int argc, char *argv[]) {
     auto const refinements = parset.get<size_t>("grid.refinements");
     grid->globalRefine(refinements);
 
+    double minDiameter = std::numeric_limits<double>::infinity();
+    double maxDiameter = 0.0;
+    for (auto it = grid->template leafbegin<0>();
+         it != grid->template leafend<0>(); ++it) {
+      auto const geometry = it->geometry();
+      auto const diam = diameter(geometry);
+      minDiameter = std::min(minDiameter, diam);
+      maxDiameter = std::max(maxDiameter, diam);
+    }
+    std::cout << "min diameter: " << minDiameter << std::endl;
+    std::cout << "max diameter: " << maxDiameter << std::endl;
+
     using GridView = Grid::LeafGridView;
     GridView const leafView = grid->leafGridView();
     size_t const leafVertexCount = leafView.size(dims);
 
+    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
+
     auto myFaces = gridConstructor.constructFaces(leafView);
 
     // Neumann boundary