06-interpolation.cc

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

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



// included standard library headers
#include <iostream>
#include <array>
#include <math.h>

// included dune-common headers
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/stringutility.hh>

// included dune-geometry headers
#include <dune/geometry/quadraturerules.hh>

// included dune-grid headers
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/uggrid.hh>

// included dune-istl headers
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>


// included dune-localfunctions headers
#include <dune/localfunctions/lagrange/pqkfactory.hh>

// included dune-fu-tutorial headers
#include <dune/fu-tutorial/referenceelementutility.hh>


// included dune-fu-tutorial headers
#include "04-gridviews.hh"
#include "05-poisson-problem.hh"



template<class GridView, class Vector, class F>
void interpolateFunction(const GridView& gridView, Vector& v, const F& f)
{
  static const int dim = GridView::dimension;
  const auto& indexSet = gridView.indexSet();
  std::size_t size = indexSet.size(dim);

  v.resize(size);
  v = 0;

  Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;

  for(const auto& e: Dune::elements(gridView))
  {
    const auto& geometry = e.geometry();
    const auto& localFiniteElement = feCache.get(e.type());
    auto localSize = localFiniteElement.localBasis().size();

    using Element = std::decay_t<decltype(e)>;
    using LocalDomain = typename Element::Geometry::LocalCoordinate;
    using Range = double;

    auto localF = [&](const LocalDomain& x) {
      return f(geometry.global(x));
    };

    std::vector<Range> localV;
    
    localFiniteElement.localInterpolation().interpolate(localF, localV);

    for(std::size_t i=0; i<localSize; ++i)
    {
      std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
      v[globalI] = localV[i];
    }
  }
}



int main(int argc, char** argv)
{
  try{
    // Maybe initialize MPI
    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);

    // Print process rank
    if(Dune::MPIHelper::isFake)
      std::cout<< "This is a sequential program." << std::endl;
    else
      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
        <<" processes!"<<std::endl;

//    using Grid = Dune::YaspGrid<2>;
//    auto gridPointer = createCubeGrid<Grid>();

//    using Grid = Dune::YaspGrid<3>;
//    auto gridPointer = createCubeGrid<Grid>();

//    using Grid = Dune::UGGrid<2>;
//    auto gridPointer = createCubeGrid<Grid>();

    using Grid = Dune::UGGrid<2>;
    auto gridPointer = createSimplexGrid<Grid>();

    Grid& grid = *gridPointer;

    grid.globalRefine(5);
    for(int i=0; i<2; ++i)
    {
      auto gridView = grid.levelGridView(grid.maxLevel());
      for(const auto& e : Dune::elements(gridView))
        if (e.geometry().center().two_norm() < .5)
          grid.mark(1, e);
      grid.adapt();
    }
    auto gridView = grid.leafGridView();
    using GridView = decltype(gridView);


    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
    using Vector = Dune::BlockVector<Dune::FieldVector<double,1>>;

    Matrix A;
    Vector rhs;
    Vector x;


    auto rhsFunction = [](auto x) {
      return (x.two_norm() < .5)*10.0;
    };
    assemblePoissonProblemPQ1(gridView, A, rhs, rhsFunction);

    x.resize(rhs.size());
    x = 0;

    auto dirichletFunction = [](auto x) {
      return sin(x[0] * 2.0*M_PI)*.1;
    };
    interpolateFunction(gridView, x, dirichletFunction);

    std::vector<bool> isBoundary;
    isBoundary.resize(rhs.size(), false);
    computeBoundaryVerticesPQ1(gridView, isBoundary);

    for(std::size_t i=0; i<rhs.size(); ++i)
    {
      if (isBoundary[i])
      {
        for(auto& entry: A[i])
          entry = 0;
        A[i][i] = 1;
        rhs[i] = x[i];
      }
    }

    Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);
    Dune::SeqILDL<Matrix,Vector,Vector> preconditioner(A,1.0);
    Dune::CGSolver<Vector> cg(op,
        preconditioner, // preconditione
        1e-4, // desired residual reduction factor
        200,   // maximum number of iterations
        2);   // verbosity of the solver

    // Object storing some statistics about the solving process
    Dune::InverseOperatorResult statistics;

    // Solve!
    cg.apply(x, rhs, statistics);

    Dune::VTKWriter<GridView> vtkWriter(gridView);
    vtkWriter.addVertexData(x, "solution");
    vtkWriter.write("06-solution");

    return 0;
  }
  catch (Dune::Exception &e){
    std::cerr << "Dune reported error: " << e << std::endl;
  }
  catch (...){
    std::cerr << "Unknown exception thrown!" << std::endl;
  }
}