1//    Discontinous Galerlin Method
2//   based on paper from
3// Riviere, Beatrice; Wheeler, Mary F.; Girault, Vivette
4// title:
5// A priori error estimates for finite element
6// methods based on discontinuous approximation spaces
7//  for elliptic problems.
8//  SIAM J. Numer. Anal. 39 (2001), no. 3, 902--931 (electronic).
9//  ---------------------------------
10//  Formulation given by Vivette Girault
11//  ------
12// Author: F. Hecht , december 2003
13// -------------------------------
14//   nonsymetric bilinear form
15//   ------------------------
16//  solve $ -\Delta u = f$ on $\Omega$ and $u= g$ on $\Gamma$
17macro dn(u) (N.x*dx(u)+N.y*dy(u) ) //  def the normal derivative
18
19mesh Th = square(10,10); // unite square
20fespace Vh(Th,P2dc);     // Discontinous P2 finite element
21fespace Xh(Th,P2);
22//  if param = 0 => Vh must be P2 otherwise we need some penalisation
23real pena=0; // a paramater to add penalisation
24varf Ans(u,v)=
25   int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)  )
26 + intalledges(Th)(//  loop on all  edge of all triangle
27       // the edge are see nTonEdge times so we / nTonEdge
28       // remark: nTonEdge =1 on border edge and =2 on internal
29       // we are in a triange th normal is the exterior normal
30       // def: jump = external - internal value; on border exter value =0
31       //      average = (external + internal value)/2, on border just internal value
32            ( jump(v)*average(dn(u)) -  jump(u)*average(dn(v))
33          + pena*jump(u)*jump(v) ) / nTonEdge
34)
35;
36func f=1;
37func g=0;
38Vh u,v;
39Xh uu,vv;
40problem A(u,v,solver="SPARSESOLVER") = Ans
41- int2d(Th)(f*v)
42- int1d(Th)(g*dn(v)  + pena*g*v)
43;
44problem A1(uu,vv,solver=CG)
45=
46 int2d(Th)(dx(uu)*dx(vv)+dy(uu)*dy(vv)) - int2d(Th)(f*vv) + on(1,2,3,4,uu=g);
47
48 A; // solve  DG
49 A1; // solve continuous
50
51plot(u,uu,cmm="Discontinue Galerkin",wait=1,value=1);
52plot(u,cmm="Discontinue Galerkin",wait=1,value=1,fill=1);
53