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

Consider the Laplace problem with homogeneous Dirichlet condition:
{ Δu = f in Ω u = 0  on  Ω

Consider a discontinuous Galerkin space V h, e.g. discontinuous P1-Lagrange space and let introduce the IP (Interior Penalty) formulation:

| find uh V h such that  Ωhuh.hv Γ {huh.n } [v ] Γ [uh ] {hv.n } +Γμ [uh ] [v ] =Ωfvv V h.

where Ω = T, Γ = T (the set of all sides of mesh elements) and, Sk, denoting the side shared by the elements Tk and T:

{v }|Sk = 1 2 ((v|Tk)|Sk + (v|T)|Sk) [v ]|Sk = ((v|Tk)|Sk (v|T)|Sk) .

For a non shared side Sk, we set

{v }|Sk = [v ]|Sk = (v|Tk)|Sk.

The operators {.} and [.] are implemented in XLiFE++ as mean(.) and jump(.) operators. The normal vector n is chosen as the outward normal from Tk on side Sk; in practice outward from the "first" parent element of the side. μ is a penalty function usually chosen to be proportional to the inverse of the measure of the side Sk. 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 :

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"); 
  BilinearForm a=intg(omega,grad(u)|grad(v))-intg(gamma,mean(grad(u)|_n)*jump(v)) 
                -intg(gamma,jump(u)*mean(grad(v)|_n))+intg(gamma,mu*jump(u)*jump(v)); 
  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 L2 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:


pict

Figure 3.1: Solution of the 2D Dirichlet problem on the unit square [0,1]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:

Ωhuh.hv Γ {huh.n } [v ] +Γ [uh ] {hv.n } +Γμ [uh ] [v ]

that is close to the IP formulation but non symmetric. For NIPG, the L2 error is weaker (about 0.00141) and the solution is less polluted at the corners.