XLiFE++ website is now hosted here with up-to-date content. The current website, with potentially deprecated content, will be removed soon.

# 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:

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.