nonlinearpoissonfem.hh

#ifndef PACXX_PROJECTSEMINAR_2019_SRC_NONLINEARPOISSONFEM_HH
#define PACXX_PROJECTSEMINAR_2019_SRC_NONLINEARPOISSONFEM_HH

#include <numeric>
#include <vector>

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

#include <dune/geometry/referenceelements.hh>

#include "lopdefaults.hh"
#include "loppattern.hh"
#include "numericaljacobian.hh"
#include "quadrature.hh"

namespace PPS {

  /** \brief a local operator for solving the nonlinear Poisson equation with
   *         conforming FEM
   *
   * \f{align*}{
   *   -\Delta u(x) + q(u(x)) &=& f(x) x\in\Omega,  \\
   *                     u(x) &=& g(x) x\in\partial\Omega_D \\
   *  -\nabla u(x) \cdot n(x) &=& j(x) x\in\partial\Omega_N \\
   * \f}
   *
   */
  template<typename Param>
  class NonlinearPoissonFEM :
    public NumericalJacobianVolume<NonlinearPoissonFEM<Param> >,
    public FullVolumePattern,
    public LocalOperatorDefaultFlags
  {
    Param& param_; // parameter functions
    int incrementorder_; // increase of integration order

  public:
    // pattern assembly flags
    enum { doPatternVolume = true };

    // residual assembly flags
    enum { doLambdaVolume = true };
    enum { doLambdaBoundary = true };
    enum { doAlphaVolume = true };

    //! constructor stores a copy of the parameter object
    NonlinearPoissonFEM (Param& param, int incrementorder=0)
      : param_(param), incrementorder_(incrementorder)
    {}

    const Param &param() const { return param_; }
    int incrementorder() const { return incrementorder_; }

    //! right hand side integral
    template<typename Element, class LocalBasis, typename R>
    void lambda_volume (const Element& elem, const LocalBasis &lb, R& r) const
    {
      // select quadrature rule
      auto geo = elem.geometry();
      const int order = incrementorder() + 2*lb.order();

      std::vector<typename LocalBasis::Traits::RangeType> phihat;

      // loop over quadrature points
      for (const auto& ip : PPS::quadratureRule(geo,order))
      {
        // evaluate basis functions
        lb.evaluateFunction(ip.position(), phihat);

        // integrate -f*phi_i
        decltype(ip.weight()) factor = ip.weight()*
          geo.integrationElement(ip.position());
        auto f=param().f(elem,ip.position());
        for (auto i : Dune::range(lb.size()))
          r[i] -= f*phihat[i]*factor;
      }
    }

    // Neumann boundary integral
    template<typename Intersection, class LocalBasis, typename R>
    void lambda_boundary (const Intersection& isect, const LocalBasis &lb,
                          R& r) const
    {
      // evaluate boundary condition type
      auto localgeo = isect.geometryInInside();
      auto facecenterlocal =
        referenceElement(localgeo).position(0,0);
      bool isdirichlet=param().b(isect, facecenterlocal);

      // skip rest if we are on Dirichlet boundary
      if (isdirichlet) return;

      // select quadrature rule
      auto globalgeo = isect.geometry();
      const int order = incrementorder() + 2*lb.order();

      std::vector<typename LocalBasis::Traits::RangeType> phihat;

      // loop over quadrature points and integrate normal flux
      for (const auto& ip : PPS::quadratureRule(globalgeo, order))
      {
        // quadrature point in local coordinates of element
        auto local = localgeo.global(ip.position());

        // evaluate shape functions (assume Galerkin method)
        lb.evaluateFunction(local, phihat);

        // integrate j
        decltype(ip.weight()) factor = ip.weight()*
          globalgeo.integrationElement(ip.position());
        auto j = param().j(isect, ip.position());
        for (auto i : Dune::range(lb.size()))
          r[i] += j*phihat[i]*factor;
      }
    }

    //! volume integral depending on test and ansatz functions
    template<typename Element, class LocalBasis, typename X, typename R>
    void alpha_volume (const Element& elem, const LocalBasis &lb, const X& x,
                       R& r) const
    {
      // types & dimension
      const int dim = Element::dimension;
      using RF = typename LocalBasis::Traits::RangeFieldType;

      // select quadrature rule
      auto geo = elem.geometry();
      const int order = incrementorder() + 2*lb.order();

      std::vector<typename LocalBasis::Traits::RangeType> phihat;
      std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
      std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;

      // loop over quadrature points
      for (const auto& ip : PPS::quadratureRule(geo, order))
      {
        // evaluate basis functions
        lb.evaluateFunction(ip.position(), phihat);

        // evaluate u
        auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));

        // evaluate gradient of shape functions
        lb.evaluateJacobian(ip.position(), gradphihat);

        // transform gradients of shape functions to real element
        const auto S = geo.jacobianInverseTransposed(ip.position());
        for (auto i : Dune::range(lb.size()))
          S.mv(gradphihat[i][0],gradphi[i][0]);

        // compute gradient of u
        Dune::FieldVector<RF,dim> gradu(0.0);
        for (auto i : Dune::range(lb.size()))
          gradu.axpy(x[i],gradphi[i][0]);

        // integrate (grad u)*grad phi_i + q(u)*phi_i
        auto factor = ip.weight()*
          geo.integrationElement(ip.position());
        auto q = param().q(u);
        for (auto i : Dune::range(lb.size()))
          r[i] += (gradu*gradphi[i][0]+
                   q*phihat[i])*factor;
      }
    }

    //! apply local jacobian of the volume term
    /**
     * \param lp The point at which to linearize
     * \param x  The point to apply the operator to
     * \param y  Where to add the result to
     */
    template<class Element, class LocalBasis, class X, class Y>
    void jacobian_apply_volume(const Element &elem, const LocalBasis &lb,
                               const X& lp, const X& x, Y& y) const
    {
      // types & dimension
      const int dim = Element::dimension;
      using RF = typename LocalBasis::Traits::RangeFieldType;

      // select quadrature rule
      auto geo = elem.geometry();
      const int order = incrementorder() + 2*lb.order();

      std::vector<typename LocalBasis::Traits::RangeType> phihat;
      std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
      std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;

      // loop over quadrature points
      for (const auto& ip : PPS::quadratureRule(geo, order))
      {
        // evaluate basis functions
        lb.evaluateFunction(ip.position(), phihat);

        // evaluate u
        auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
        // w is the function corresponding to coefficients lp
        auto w = std::inner_product(begin(lp), end(lp), begin(phihat), RF(0));

        // evaluate gradient of shape functions
        lb.evaluateJacobian(ip.position(), gradphihat);

        // transform gradients of shape functions to real element
        const auto S = geo.jacobianInverseTransposed(ip.position());
        for (auto i : Dune::range(lb.size()))
          S.mv(gradphihat[i][0],gradphi[i][0]);

        // compute gradient of u
        Dune::FieldVector<RF,dim> gradu(0.0);
        for (auto i : Dune::range(lb.size()))
          gradu.axpy(x[i],gradphi[i][0]);

        // integrate (grad u)*grad phi_i + q(u)*phi_i
        auto factor = ip.weight()*
          geo.integrationElement(ip.position());
        auto qprime_u = param().qprime(w)*u;
        for (auto i : Dune::range(lb.size()))
          y[i] += (gradu*gradphi[i][0]+
                   qprime_u*phihat[i])*factor;
      }
    }

  };

}

#endif // PACXX_PROJECTSEMINAR_2019_SRC_NONLINEARPOISSONFEM_HH