I don’t know if you checked the L2 from the example DC but it is pretty close to

the paper although their method appears slightly better its not much.

I can’t get their method to work yet but you may want to look at this for comparison.

See if this helps. I think though the error along x=y may be informative. I’m doing this

to test mjmdatascope but there may be a similar way with existing FF code.

Obvoously you can just delete that or I can put everything up on github if there is interest.

```
load "msh3"
load "medit"
load "mjm_dscope_freefem"
// mjm - I added the formulation from this,
// and it is very close but slightly worse than the paper
// https://www.sciencedirect.com/science/article/pii/S0377042714002131#fd000050
// although I still can't get that to work. ZZ
// https://www.um.es/freefem/ff++/pmwiki.php?n=Main.DiscontinuousGalerkin
// Discontinous Galerlin Method
// based on paper from
// Riviere, Beatrice; Wheeler, Mary F.; Girault, Vivette
// title:
// A priori error estimates for finite element
// methods based on discontinuous approximation spaces
// for elliptic problems.
// SIAM J. Numer. Anal. 39 (2001), no. 3, 902--931 (electronic).
// ---------------------------------
// Formulation given by Vivette Girault
// ------
// Author: F. Hecht , december 2003
// -------------------------------
macro dn(u) (N.x*dx(u)+N.y*dy(u) ) // def the normal derivative
int nref = 6;
real [int] L2errorm ( nref ); // initialize the L2 error d array
for ( int nn = 0; nn < nref; nn++) {
int Nn = 2^(nn + 1) ; // space discretization
// "2" means 2 squares or a side of h but hmax has sqrt(2)
real h = 1.0/Nn;
mesh Th = square (Nn,Nn); // mesh generation of a square
//medit(""+Nn,Th);
//fespace Vh(Th,P2); // piecewise continuous quadratic finite element.
//fespace Vh(Th,P1); // Discontinous P2 finite element
fespace Vh(Th,P2dc); // Discontinous P2 finite element
Vh uh,vh;
Vh uex = sin(pi*x) * sin (pi*y); // exact solution
Vh f = 2.0 * sin(pi*x) *(pi^2) * sin (pi*y); // corresponding RHS
//
varf a(uh,vh, solver=sparsesolver,tgv=-1) = int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh) )
// this seems to be important? for the continuous case anyway
+on(0,1,2,3,4,uh=0)
+intalledges(Th)( ( jump(vh)*mean(dn(uh)) - jump(uh)*mean(dn(vh))
// + pena*jump(u)*jump(v)
) / nTonEdge)
// + intalledges(Th)(sqrt(1)*(1.0/Nn)^-1 * jump(uh)*jump(vh) ); // bilinear form
// + intalledges(Th)((lenEdge)^-1 * jump(uh)*jump(vh) ); // bilinear form
;
matrix A = a(Vh, Vh,tgv=-1); // build the matrix
//set(A,solver=CG);
set(A,solver=sparsesolver);
varf L(uh, vh,tgv=-1) = int2d (Th)(f*vh)+ on(0,1,2,3,4,uh=0) ; // linear form
Vh F; F[] = L(0, Vh,tgv=-1); // build the right hand side vector
uh[] = A^-1 * F [];
int na=uh[].n;
//real s=0;
//for(int i=0; i<na; ++i) s=s+(uh[][i]-uex[][i])^2;
//L2errorm [nn] = sqrt(s)*h*h; // sqrt ( int2d (Th)((uh - uex )^2));
L2errorm [nn] = sqrt ( int2d (Th)((uh - uex )^2));
meshL Th1=Sline(1024,[x]);
fespace Mh(Th1,P1);
Mh asol=uh(x,x);
Mh anex=uex(x,x);
Mh diff=anex-asol; // =uex(x,x);
mjmdscopeL(Th1,"names asol anex diff "+nn,asol,anex,diff); // warning do not forget ()
}
//
for ( int n = 0; n < nref; n ++){
int Nn = 2^(n + 1);
cout << " N " << " = " << Nn << endl;
cout << " L2errorm " << n << " = " << L2errorm [n] << endl;}
for ( int n = 1; n < nref; n ++) {
int Nn = 2^(n + 1);
cout << " N " << " = " << Nn << endl;
cout <<" convergence rate d = " << log ( L2errorm [n - 1]/ L2errorm [n])/log (2.) << endl; }
```

MWGFEM_1.edp (3.1 KB)