Bassi-rebay mixed formulation

General Discussion
1

Context

I want to solve the heat equation with a (Discontinuous Galerkin) Bassi-Rebay formulation. My toy problem is -\Delta u = -4\pi² sin(2\pi x) \quad \forall x \in [0,1] with homogeneous Dirichlet conditions. I am implementing a 1D solver (with Legendre polynomials), but I have an conditioning issue is my square matrix of the linear system. So I first want to try a FreeFem implementation before continuing.

For remainder the Bassi-Rebay formulation is:

\begin{align}a(q_h, r_h) + b(u_h, r_h) &= F(r_h)\\ -b(v_h, q_h)&= G(v_h) \end{align}

where:

\begin{align}a(q_h, r_h) &= \sum_K \int_K q_h\cdot r_h\\ b(u_h,r_h) &= -\sum_K \int_K \nabla u_h\cdot r_h + \sum_{\Gamma^{ID}}\int_\Gamma \braket{r_h}\cdot n[u_h]\\ F(r_h) &= \sum_{\Gamma^D} \int_\Gamma u_Dr_h\cdot n\\ G(v_h) &= \int_\Omega fv_h + \sum_{\Gamma^N} \int_\Gamma g_N v_h \end{align}

With my current FreeFem script, I get some weird results. Here is my input file:

mesh th = square(10, 1, [2*x - 1, y]);

fespace Xh(th, P1dc);
Xh u, v, q, r;

func f = -4 * pi^2 * sin(2 * pi * x);

problem BassiRebay([q, u], [r, v]) =
    int2d(th)( q * r )

    - int2d(th)( dx(u) * r )
    + int2d(th)( q * dx(v) )

    + intalledges(th)( (abs(N.x) > 0.5) * mean(r) * jump(u) )
    - intalledges(th)( (abs(N.x) > 0.5) * mean(q) * jump(v) )
    + on(4, 2, u = 0)
    - int2d(th)( f * v );

BassiRebay;

plot(u, fill=1, value=1, wait=1, cmm="Bassi-Rebay with jump/mean operators");
2

Your code solves the problem as a 2d problem. In order to have a well-posed problem (non singular matrix) you have to include the y derivatives, with for example Neumann boundary condition at the bottom and top.
Moreover, there was a sign error for the jump terms. This is probably related to the FreeFem convention that the jump is “external minus internal” whereas it can be the opposite in several papers.

mesh th = square(10, 1);

fespace Xh(th, P1dc);
Xh u, v, qx,qy, rx,ry;

func f = 4 * pi^2 * sin(2 * pi * x);

problem BassiRebay([u, qx, qy], [v, rx, ry]) =
    int2d(th)( qx * rx + qy * ry )

    - int2d(th)( dx(u) * rx + dy(u) * ry )
    + int2d(th)( qx * dx(v) + qy * dy(v) )

    - intalledges(th)( (mean(rx)*N.x+mean(ry)*N.y) * jump(u)*(nTonEdge-1)/2. )
    + intalledges(th)( (mean(qx)*N.x+mean(qy)*N.y) * jump(v)*(nTonEdge-1)/2. )

    + int1d(th,2,4)((rx*N.x+ry*N.y) * u)
    - int1d(th,2,4)((qx*N.x+qy*N.y) * v)

    - int2d(th)( f * v );

BassiRebay;

real error=sqrt(int2d(th)((u-sin(2.*pi*x))^2));
cout << "L2 error = " << error << endl;

plot(u, fill=1, value=1, wait=1, dim=3, cmm="Bassi-Rebay with jump/mean operators");