/*
@book{Q12,
  address =       {Italy},
  author =        {A. Quarteroni},
  edition =       {5th},
  publisher =     {Springer-Verlag},
  title =         {Modellistica numerica per problemi differenziali},
  year =          {2012},
  isbn =          {978-88-470-2747-3},
}

page 355
*/
// -mu \nabla^2 u + b * \nabla u = f \Omega
// u = g                             \Gamma
real mu = 1e-1,b1 = 1.0, b2 = 1.0;
macro LS(u) (-mu*(dxx(u)+dyy(u))) // symmetric
macro LSS(u) (b1*dx(u)+b2*dy(u)) // skew-symmetric
macro L(u) LS(u)+LSS(u) // linear operator
macro dn(u) (N.x*dx(u)+N.y*dy(u)) //  def the normal derivative
func uexact = x+y*(1-x)+(exp(-1.0/mu)-exp(-(1-x)*(1-y)/mu))/(1-exp(-1.0/mu));
func g = x+y*(1-x)+(exp(-1.0/mu)-exp(-(1-x)*(1-y)/mu))/(1-exp(-1.0/mu));
func f = exp(-(1-x)*(1-y)/mu)/(mu*(1-exp(-1.0/mu)))*
  ((x-1)^2+(y-1)^2+(x-1)+(y-1))-x-y+2;
int n = 16;
mesh Th = square(n,n); // unit square
fespace Xh(Th,P1);
real delta = 0.5; // (0 Galerkin)
real sigma = 0.; // (1 GLS, 0 SUPG, -1 DW)
Xh uh,vh;
solve ADG(uh,vh) = int2d(Th)(mu*dx(uh)*dx(vh)+mu*dy(uh)*dy(vh))+
  int2d(Th)((b1*dx(uh)+b2*dy(uh))*vh)+
  int2d(Th)(LSS(uh)*delta*hTriangle/sqrt(b1^2+b2^2)*LSS(vh))+
  int2d(Th)(LS(uh)*delta*hTriangle/sqrt(b1^2+b2^2)*LSS(vh))-
  int2d(Th)(f*delta*hTriangle/sqrt(b1^2+b2^2)*LSS(vh))+
  int2d(Th)(LSS(uh)*delta*hTriangle/sqrt(b1^2+b2^2)*sigma*LS(vh))+
  int2d(Th)(LS(uh)*delta*hTriangle/sqrt(b1^2+b2^2)*sigma*LS(vh))-
  int2d(Th)(f*delta*hTriangle/sqrt(b1^2+b2^2)*sigma*LS(vh))-
  int2d(Th)(f*vh)+on(1,2,3,4,uh=g);
plot(uh,cmm="(Generalized) Galerkin",wait=1,value=1,fill=1);
cout << "Generalized Galerkin err: " <<
  sqrt(int2d(Th)((uh-uexact)^2)) << endl;
fespace Wdelta(Th,P1dc); int r=1;    // Discontinous elements
Wdelta udelta,vdelta,gdelta;
gdelta = g;
real gamma=10.0*r^2, tau=1.0;
solve ADDG(udelta,vdelta)=
  int2d(Th)(mu*dx(udelta)*dx(vdelta)+mu*dy(udelta)*dy(vdelta))
  +intalledges(Th)(mu*jump(vdelta)*mean(dn(udelta))/nTonEdge)
  +intalledges(Th)(tau*mu*jump(udelta)*mean(dn(vdelta))/nTonEdge)
  +int1d(Th)(tau*mu*gdelta*dn(vdelta))
  +intalledges(Th)(gamma/lenEdge*jump(udelta)*jump(vdelta)/nTonEdge)
  -int1d(Th)(gamma/lenEdge*gdelta*vdelta)
  -int2d(Th)(udelta*(b1*dx(vdelta)+b2*dy(vdelta)))
  +intalledges(Th)((abs(b1*N.x+b2*N.y)*jump(udelta)/2-
                    (b1*N.x+b2*N.y)*mean(udelta)/(3.0-nTonEdge))
                   *jump(vdelta)/nTonEdge)
  +int1d(Th)(min(b1*N.x+b2*N.y,0.0)*gdelta*vdelta)
  -int2d(Th)(f*vdelta);
plot(udelta,cmm="Discontinuous Galerkin",wait=1,value=1,fill=1);
cout << "Discontinuous Garlerkin err: " <<
  sqrt(int2d(Th)((udelta-uexact)^2)) << endl;