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::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);
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()))
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)
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
//! 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