Newer
Older
#ifndef PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_HH
#define PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_HH
#include <dune/common/fvector.hh>
#include <dune/common/rangeutilities.hh>
namespace PPS {
/** \brief a local operator for solving the nonlinear Poisson equation with
*
* \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 NumericalJacobianBoundary<NitscheNonlinearPoissonFEM<Param> >
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::RangeField;
// 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::array<typename LocalBasis::Jacobian, lb.size()> gradphi;
// 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
auto phihat = lb.evaluateFunction(local);
auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
// evaluate gradient of shape functions
auto gradphihat = lb.evaluateJacobian(local);
// 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()))
auto n = isect.unitOuterNormal(ip.position());
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()))
(-(gradu*n)*phihat[i]
-(u-g)*(gradphi[i][0]*n)
//! 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::RangeField;
// 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::array<typename LocalBasis::Jacobian, lb.size()> gradphi;
// 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
auto phihat = lb.evaluateFunction(local);
// evaluate u
auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
// evaluate gradient of shape functions
auto gradphihat = lb.evaluateJacobian(local);
// 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