diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc index 0dfdc6f0d0de350db42e9704186387fd442f953e..cf8318988177a3b2af04a1b64eb0f358b77701c5 100644 --- a/src/sand-wedge.cc +++ b/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