# Example 7/18 - Discontinuous Galerkin method for 2D Dirichlet problem

Consider the Laplace problem with homogeneous Dirichlet condition:

Consider a discontinuous Galerkin space ${V}_{h}$, e.g. discontinuous ${P}^{1}$-Lagrange space and let introduce the IP (Interior Penalty) formulation:

where $\mathrm{\Omega }=\bigcup {T}_{\ell }$, $\Gamma =\bigcup \partial {T}_{\ell }$ (the set of all sides of mesh elements) and, ${S}_{k,\ell }$ denoting the side shared by the elements ${T}_{k}$ and ${T}_{\ell }$:

$\left\{v{\right\}}_{|{S}_{k\ell }}=\frac{1}{2}\left({\left({v}_{|{T}_{k}}\right)}_{|{S}_{k\ell }}+{\left({v}_{|{T}_{\ell }}\right)}_{|{S}_{k\ell }}\right)\phantom{\rule{1em}{0ex}}\left[v{\right]}_{|{S}_{k\ell }}=\left({\left({v}_{|{T}_{k}}\right)}_{|{S}_{k\ell }}-{\left({v}_{|{T}_{\ell }}\right)}_{|{S}_{k\ell }}\right).$

For a non shared side ${S}_{k}$, we set

$\left\{v{\right\}}_{|{S}_{k}}=\left[v{\right]}_{|{S}_{k}}={\left({v}_{|{T}_{k}}\right)}_{|{S}_{k}}.$

The operators $\left\{.\right\}$ and $\left[.\right]$ are implemented in XLiFE++ as mean(.) and jump(.) operators. The normal vector $n$ is chosen as the outward normal from ${T}_{k}$ on side ${S}_{k\ell }$; in practice outward from the "first" parent element of the side. $\mu$ is a penalty function usually chosen to be proportional to the inverse of the measure of the side ${S}_{k\ell }$. Note that the Dirichlet boundary condition is a natural condition in this formulation.

To deal with a such formulation, the following objects have to be constructed from a geometrical domain, say omega :

• a FE space specifying $L2$ conformity (discontinuous space) defined on omega

• the geometrical domain Gamma of sides of omega using the function sides(omega)

• the bilinear form related to IP formulation and involving mean(.) and jump(.) operators

The XLiFE++ implementation of this problem using discontinuous P1-Lagrange is the following:

#include "xlife++.h"
using namespace xlifepp;

Real uex(const Point& p, Parameters& pars=defaultParameters)
{return sin(pi_*p(1))*sin(pi_*p(2));}
Real f(const Point& p, Parameters& pars=defaultParameters)
{return 2*pi_*pi_*sin(pi_*p(1))*sin(pi_*p(2));}
Real fmu(const Point& p, Parameters& pars=defaultParameters)
{GeomElement* gelt=getElementP();
if(gelt!=0) return 1/gelt->measure();
return 0.;
}

int main(int argc, char** argv)
{
init(argc, argv, _lang=en); // mandatory initialization of xlifepp
// create mesh and geometrical domains
Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=31,_domain_name="Omega");
Mesh mR(R,_triangle,1,_structured,_alternateSplit);
Domain omega=mR.domain("Omega");
Domain gamma=sides(omega);  // create domain of sides
// create discontinuous P1Lagrange space
Interpolation in(Lagrange,standard,1,L2);
Space V(omega,in,"V");
Unknown u(V,"u");TestFunction v(u,"v");
// create bilinear form and TermMatrix
Function mu(fmu); mu.require("element");
TermMatrix A(a,"A");
// create the right hand side an solve the linear system
TermVector F(intg(omega,f*v),"F");
TermVector U=directSolve(A,F);
saveToFile("U",U,_vtu);
// compute L2 error
TermMatrix M(intg(omega,u*v),"M");
TermVector Uex(u,omega,uex);
TermVector E=U-Uex;
theCout<<" IP : |Uuex|L2 = "<<sqrt(M*E|E)<<eol;
return 0;
}

The ${L}^{2}$ error is about $0.00317$ for 5400 dofs. Using Paraview with the Warp by scalar filter that produces elevation, the approximated field $u$ looks like: Figure 3.1: Solution of the 2D Dirichlet problem on the unit square ${\left[0,1\right]}^{2}$ with discontinuous P1-Lagrange elements, IP formulation (left) and NIPG formulation (right)

It can be noticed that the field is discontinuous and major errors is located on the corners. NIPG formulation is an other penalty method corresponding to the bilinear form:

${\int }_{\mathrm{\Omega }}{\nabla }_{h}{u}_{h}.{\nabla }_{h}v-{\int }_{\Gamma }\left\{{\nabla }_{h}{u}_{h}.n\right\}\phantom{\rule{0.17em}{0ex}}\left[v\right]+{\int }_{\Gamma }\left[{u}_{h}\right]\phantom{\rule{0.17em}{0ex}}\left\{{\nabla }_{h}v.n\right\}+{\int }_{\Gamma }\mu \left[{u}_{h}\right]\phantom{\rule{0.17em}{0ex}}\left[v\right]$

that is close to the IP formulation but non symmetric. For NIPG, the ${L}^{2}$ error is weaker (about $0.00141$) and the solution is less polluted at the corners.