nitschenonlinearpoissonfem.hh

#ifndef PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_HH
#define PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_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 "nonlinearpoissonfem.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 NitscheNonlinearPoissonFEM :
    public NonlinearPoissonFEM<Param>,
    public NumericalJacobianBoundary<NitscheNonlinearPoissonFEM<Param> >
  {
    using Real = typename Param::value_type;
    Real stab_;

  public:
    // residual assembly flags
    enum { doAlphaBoundary = true };

    //! constructor stores a copy of the parameter object
    NitscheNonlinearPoissonFEM (Param& param, Real stab,
                                int incrementorder=0) :
      NonlinearPoissonFEM<Param>{param, incrementorder},
      stab_{stab}
    {}

    //! volume integral depending on test and ansatz functions
    template<typename Intersection, class LocalBasis, typename X, typename R>
    void alpha_boundary (const Intersection& isect, const LocalBasis &lb,
                         const X& x, R& r) const
    {
      // evaluate boundary condition type
      auto localgeo = isect.geometryInInside();
      auto facecenterlocal =
        referenceElement(localgeo).position(0,0);
      bool isdirichlet=this->param().b(isect, facecenterlocal);

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

      // types & dimension
      const int dim = Intersection::Entity::dimension;
      using RF = typename LocalBasis::Traits::RangeFieldType;

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

      // geometry of inside cell
      auto geo_inside = isect.inside().geometry();

      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))
      {
        // quadrature point in local coordinates of element
        auto local = localgeo.global(ip.position());

        // evaluate basis functions
        lb.evaluateFunction(local, phihat);

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

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

        // transform gradients of shape functions to real element
        const auto S = geo_inside.jacobianInverseTransposed(local);
        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]);

        // get unit outer normal vector
        auto n = isect.unitOuterNormal(ip.position());
        // get dirichlet boundary condition value
        auto g = this->param().g(isect, ip.position());

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

    //! volume integral depending on test and ansatz functions
    template<typename Intersection, class LocalBasis, typename X, typename Y>
    void jacobian_apply_boundary(const Intersection& isect,
                                 const LocalBasis &lb,
                                 const X& lp, const X& x, Y& y) const
    {
      // evaluate boundary condition type
      auto localgeo = isect.geometryInInside();
      auto facecenterlocal =
        referenceElement(localgeo).position(0,0);
      bool isdirichlet=this->param().b(isect, facecenterlocal);

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

      // types & dimension
      const int dim = Intersection::Entity::dimension;
      using RF = typename LocalBasis::Traits::RangeFieldType;

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

      // geometry of inside cell
      auto geo_inside = isect.inside().geometry();

      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))
      {
        // quadrature point in local coordinates of element
        auto local = localgeo.global(ip.position());

        // evaluate basis functions
        lb.evaluateFunction(local, phihat);

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

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

        // transform gradients of shape functions to real element
        const auto S = geo_inside.jacobianInverseTransposed(local);
        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]);

        // get unit outer normal vector
        auto n = isect.unitOuterNormal(ip.position());

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

}

#endif // PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_HH