driver.hh

#ifndef PACXX_PROJECTSEMINAR_2019_SRC_DRIVER_HH
#define PACXX_PROJECTSEMINAR_2019_SRC_DRIVER_HH

#include <cstddef>
#include <iostream>
#include <memory>
#include <set>
#include <utility>

#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>

#include <dune/istl/bvector.hh>

#include <dune/geometry/refinement.hh>

#include <dune/grid/io/file/vtk/common.hh>
#include <dune/grid/io/file/vtk/function.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>

#include "constraints.hh"
#include "gridoperator.hh"
#include "interpolate.hh"
#include "newton.hh"
#include "problem.hh"
#include "vtkfunction.hh"

#ifdef NITSCHE
#include "nitschenonlinearpoissonfem.hh"
#else
#include "nonlinearpoissonfem.hh"
#endif

namespace PPS {

  template<typename GV>
  void driver (const GV& gv, Dune::ParameterTree& ptree)
  {
    // dimension and important types
    typedef double RF;                   // type for computations

    // make PDE parameter class
    RF eta = ptree.get("problem.eta",(RF)1.0);
    using Problem = NonlinearPoissonProblem<RF>;
    Problem problem(eta);
    auto glambda = [&](const auto& e, const auto& x)
      {return problem.g(e,x);};

    // Assemble constraints
    //== Exercise 2 {
    std::set<std::size_t> cc;
#ifndef NITSCHE
    auto isDirichlet = [&](const auto& intersection) {
      auto testpoint = referenceElement(intersection.geometry()).position(0,0);
      return problem.b(intersection,testpoint);
    };
    PPS::p1Constraints(isDirichlet,gv, cc); // assemble constraints
    std::cout << "constrained dofs=" << cc.size() << " of "
              << gv.indexSet().size(gv.dimension) << std::endl;
#endif
    //== }

    // A coefficient vector
    using Z = Dune::BlockVector<Dune::FieldVector<RF, 1> >;
    Z z(gv.indexSet().size(gv.dimension)); // initial value

    // Fill the coefficient vector
    PPS::interpolate(gv,glambda,z);

    // Make a local operator
    //== Exercise 2 {
#ifdef NITSCHE
    RF stab = ptree.get("fem.stab",(RF)1);
    typedef NitscheNonlinearPoissonFEM<Problem> LOP;
    LOP lop(problem,stab);
#else
    typedef NonlinearPoissonFEM<Problem> LOP;
    LOP lop(problem);
#endif
    //== }

    // Make a global operator
    typedef PPS::GridOperator<
      GV,
      LOP,      /* local operator */
      RF  /* range/jacobian field type*/
      > GO;
    //== Exercise 2 {
    GO go(gv, lop, std::move(cc));
    //== }

    // Select a linear solver backend
    auto ls = makeSolver<Z>(ptree.sub("solver"), gv.dimension);

    // set up nonlinear solver
    PPS::Newton<GO,PPS::SolverBackend<Z>,Z> newton(go,z,*ls);
    newton.setReassembleThreshold(0.0); // always reassemble J
    newton.setVerbosityLevel(3);        // be verbose
    newton.setReduction(1e-10);         // total reduction
    newton.setMinLinearReduction(1e-4); // min. red. in lin. solve
    newton.setMaxIterations(25);        // limit number of its
    newton.setLineSearchMaxIterations(10); // limit line search
    newton.setParameters(ptree.sub("newton"));

    // solve nonlinear problem
    newton.apply();

    // Write VTK output file
    int subsampling = ptree.get("output.subsampling",(int)1);
    Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,Dune::refinementIntervals(subsampling));
    using Function = Dune::P1VTKFunction<GV, Z>;
    vtkwriter.addVertexData(std::make_shared<Function>(gv, z,"fesol"));
    vtkwriter.addVertexData(PPS::makeScalarVTKGridFunction<GV>(glambda, "exact"));
    vtkwriter.write(ptree.get("output.filename","output"),
                    Dune::VTK::appendedraw);
  }

}

#endif // PACXX_PROJECTSEMINAR_2019_SRC_DRIVER_HH