convergence.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <iostream>

#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/test/testsuite.hh>
#include <dune/curvedsurfacegrid/curvedsurfacegrid.hh>
#include <dune/curvedsurfacegrid/geometries/ellipsoid.hh>
#include <dune/curvedsurfacegrid/geometries/sphere.hh>
#include <dune/curvedsurfacegrid/geometries/torus.hh>
#include <dune/foamgrid/foamgrid.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/localfunctions/lagrange/pk.hh>
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>

// 1 .. Sphere
// 2 .. Ellipsoid
// 3 .. Torus
#ifndef GEOMETRY_TYPE
#define GEOMETRY_TYPE 1
#endif

const int order = 3;
const int quad_order = order+6;

#if GEOMETRY_TYPE < 3
const int num_levels = 4;
#else
const int num_levels = 3;
#endif

// Test Hausdorf-distance
template <class Grid>
typename Grid::ctype inf_error (const Grid& grid)
{
  std::cout << "  inf_error..." << std::endl;
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  typename Grid::ctype dist = 0;
  auto projection = grid.gridFunction().f();
  for (const auto& element : elements(grid.leafGridView()))
  {
    auto geometry = element.geometry();

    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
      auto x = geometry.global(qp.position());
      auto y = projection(x);
      dist = std::max(dist, (x - y).two_norm());
    }
  }

  return dist;
}

// Test integrated Hausdorf-distance
template <class Grid>
typename Grid::ctype L2_error (const Grid& grid)
{
  std::cout << "  L2_error..." << std::endl;
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  typename Grid::ctype dist = 0;
  auto projection = grid.gridFunction().f();
  for (const auto& element : elements(grid.leafGridView()))
  {
    auto geometry = element.geometry();

    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
      auto x = geometry.global(qp.position());
      auto y = projection(x);
      dist += (x - y).two_norm2() * qp.weight() * geometry.integrationElement(qp.position());
    }
  }

  return std::sqrt(dist);
}

// Test integrated error in normal vectors
template <class Grid>
typename Grid::ctype normal_error (const Grid& grid)
{
  std::cout << "  normal_error..." << std::endl;
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  typename Grid::ctype dist = 0;

  auto projection = grid.gridFunction().f();
  for (const auto& element : elements(grid.leafGridView()))
  {
    auto geometry = element.geometry();
    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
      auto x = geometry.impl().normal(qp.position());
      auto y = projection.normal(geometry.global(qp.position()));
      dist += (x - y).two_norm2() * qp.weight() * geometry.integrationElement(qp.position());
    }
  }

  return std::sqrt(dist);
}

// Test integrated error in mean curvature
template <class Grid>
typename Grid::ctype curvature_error (const Grid& grid)
{
  std::cout << "  curvature_error..." << std::endl;
  using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  using FiniteElementType = Dune::PkLocalFiniteElement<typename Grid::ctype, typename Grid::ctype, 2, order+4>;
  FiniteElementType fe;

  using LBTraits = typename FiniteElementType::Traits::LocalBasisType::Traits;
  std::vector<typename LBTraits::JacobianType> shapeGradients;
  std::vector<GlobalCoordinate> gradients;
  std::vector<GlobalCoordinate> normals;

  auto projection = grid.gridFunction().f();

  typename Grid::ctype dist = 0;
  for (const auto& element : elements(grid.leafGridView()))
  {
    auto geometry = element.geometry();
    const auto& localBasis = fe.localBasis();

    { // get coefficients of normal vectors
      const auto& localInterpolation = fe.localInterpolation();

      auto n = [&](const auto& local) -> GlobalCoordinate
      {
        auto x = geometry.impl().normal(local);
        return x;
      };
      localInterpolation.interpolate(n, normals);
    }

    shapeGradients.resize(localBasis.size());
    gradients.resize(localBasis.size());

    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
      auto jInvT = geometry.jacobianInverseTransposed(qp.position());
      localBasis.evaluateJacobian(qp.position(), shapeGradients);

      for (std::size_t i = 0; i < gradients.size(); ++i)
        jInvT.mv(shapeGradients[i][0], gradients[i]);

      typename Grid::ctype H = 0;
      for (std::size_t i = 0; i < gradients.size(); ++i)
        H += normals[i].dot(gradients[i]);

      auto H_exact = projection.mean_curvature(geometry.global(qp.position()));
      dist += std::abs(H/2 - H_exact) * qp.weight() * geometry.integrationElement(qp.position());
    }
  }

  return dist;
}

// longest edge in the grid
template <class Grid>
typename Grid::ctype edge_length (const Grid& grid)
{
  typename Grid::ctype h = 0;
  for (const auto& e : edges(grid.hostGrid().leafGridView()))
    h = std::max(h, e.geometry().volume());

  return h;
}


int main (int argc, char** argv)
{
  using namespace Dune;
  MPIHelper::instance(argc, argv);

  auto grid_order = std::integral_constant<int,order>{};

  using HostGrid = FoamGrid<2,3>;
#if GEOMETRY_TYPE == 1
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "sphere.msh");
  auto gridFct = sphereGridFunction<HostGrid>(1.0);
#elif GEOMETRY_TYPE == 2
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "ellipsoid3.msh");
  auto gridFct = ellipsoidGridFunction<HostGrid>(1.0, 0.75, 1.25);
#elif GEOMETRY_TYPE == 3
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "torus.msh");
  auto gridFct = torusGridFunction<HostGrid>(2.0, 1.0);
#endif

  CurvedSurfaceGrid grid(*hostGrid, gridFct, grid_order);
  using Grid = decltype(grid);
  std::vector<typename Grid::ctype> inf_errors, L2_errors, normal_errors, curvature_errors, edge_lengths;
  for (int i = 0; i < num_levels; ++i) {
    if (i > 0)
      grid.globalRefine(1);

    using DC = LagrangeDataCollector<typename Grid::LeafGridView, order>;
    VtkUnstructuredGridWriter<typename Grid::LeafGridView, DC> writer(grid.leafGridView());
    writer.write("level_" + std::to_string(i) + ".vtu");

    std::cout << "level = " << i << std::endl;
    inf_errors.push_back( inf_error(grid) );
    L2_errors.push_back( L2_error(grid) );
    normal_errors.push_back( normal_error(grid) );
    curvature_errors.push_back( curvature_error(grid) );
    edge_lengths.push_back( edge_length(grid) );
  }

  using std::log;
  using std::abs;

  // calculate the experimental order of convergence
  std::vector<typename Grid::ctype> eocInf(num_levels, 0), eocL2(num_levels, 0), eocNormal(num_levels, 0), eocCurvature(num_levels, 0);
  for (int i = 1; i < num_levels; ++i) {
    eocInf[i]       = log(inf_errors[i]/inf_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocL2[i]        = log(L2_errors[i]/L2_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocNormal[i]    = log(normal_errors[i]/normal_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocCurvature[i] = log(curvature_errors[i]/curvature_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
  }

  TestSuite test;
  test.check(abs(eocInf.back()    - (order+1)) < 0.5, "|d|_L^inf != order+1");
  test.check(abs(eocL2.back()     - (order+1)) < 0.5, "|d|_L^2 != order+1");
  test.check(abs(eocNormal.back() - (order+0)) < 0.5, "|n - nh|_L^2 != order");
  test.check(abs(eocCurvature.back()   - (order-1)) < 0.5, "|H - Hh|_L^2 != order-1");


  auto print_line = [](auto i, const auto& data, auto accept) {
    std::cout << i;
    for (const auto& d : data) {
      std::cout << '\t' << "| ";
      if (accept(d))
        std::cout << d;
      else
        std::cout << '\t';
    }
    std::cout << std::endl;
  };

  auto print_break = [] {
    std::cout.width(8 + 8*16);
    std::cout.fill('-');
    std::cout << '-' << std::endl;
    std::cout.fill(' ');
  };

  std::cout.setf(std::ios::scientific);
  print_line("level", std::vector<std::string>{"err_inf", "eoc_inf", "err_L2", "eoc_L2", "err_norm", "eoc_norm", "err_curv", "eoc_curv"}, [](auto&&) { return true; });
  print_break();

  for (int i = 0; i < num_levels; ++i)
    print_line(i, std::vector<typename Grid::ctype>{inf_errors[i], eocInf[i], L2_errors[i], eocL2[i], normal_errors[i], eocNormal[i], curvature_errors[i], eocCurvature[i]},
      [](auto const& v) { return v != 0.0; });
  print_break();


  return test.exit();
}