// -\Delta u = f \Omega=[0,1]^2 // u = g \partial\Omega // macro dn(u) (N.x*dx(u)+N.y*dy(u) ) // def the normal derivative func f=1; //func g=1.0*(x==0)*(y<1)+1.0*(y==1); real alpha = 10; func g = (exp(alpha*(x+y))-1)/(exp(2*alpha)-1); mesh Th = square(20,20); // unit square fespace Wdelta(Th,P1dc); // Discontinous P1 finite element cout << "DG dof " << Wdelta.ndof << endl; Wdelta udelta,vdelta,gdelta; gdelta = g; real gamma=1, tau=-1; solve PoissonDG(udelta,vdelta)= int2d(Th)(dx(udelta)*dx(vdelta)+dy(udelta)*dy(vdelta)) +intalledges(Th)((mean(dn(udelta))*jump(vdelta) +tau*jump(udelta)*mean(dn(vdelta)) +gamma/lenEdge*jump(udelta)*jump(vdelta))/nTonEdge) +int1d(Th)(tau*gdelta*dn(vdelta)-gamma/lenEdge*gdelta*vdelta)-int2d(Th)(f*vdelta); plot(udelta,cmm="DG-N",wait=1,value=1,fill=1); fespace Vh(Th,P1); cout << "Galerkin dof " << Vh.ndof << endl; Vh uh,vh; solve Poisson(uh,vh)=int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh)) -int2d(Th)(f*vh)+on(1,2,3,4,uh=g); plot(uh,cmm="Galerkin",wait=1,value=1,fill=1);