# Bassi-rebay mixed formulation
**URL:** https://community.freefem.org/t/bassi-rebay-mixed-formulation/4237
**Category:** General Discussion
**Created:** [March 8, 2026, 8:45pm UTC](https://community.freefem.org/t/bassi-rebay-mixed-formulation/4237 "2026-03-08T20:45:52Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@MaximeLee](https://community.freefem.org/u/MaximeLee)
#### Post date: [March 8, 2026, 8:45pm UTC](https://community.freefem.org/t/bassi-rebay-mixed-formulation/4237/1 "2026-03-08T20:45:52Z")
# 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:
```auto
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");
```
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [March 9, 2026, 11:35am UTC](https://community.freefem.org/t/bassi-rebay-mixed-formulation/4237/2 "2026-03-09T11:35:07Z")
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.
```auto
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");
```