Divergence-free condition on boundary

Hello everyone, I have a question about boundary conditions. I want to solve the incompressible MHD equations, as follows:

For the magnetic field, the boundary conditions I chose are:
\mathbf{B}\times\mathbf{n}|_{\partial \Omega}=0,\quad \nabla\cdot \mathbf{B}|_{\partial \Omega}=0

Regarding the boundary conditions \mathbf{B}\times\mathbf{n}|_{\partial \Omega}=0, I know the correct approach is:

+ int1d(Th,1,2,3,4)(1e10*crossuB(N.x,N.y,Bx,By)*crossuB(N.x,N.y,cx,cy))

with

macro crossuB(ux,uy,Bx,By) ((ux)*(By) - (uy)*(Bx))

However, I applied a divergence-free condition to the boundary \nabla\cdot \mathbf{B}|_{\partial \Omega}=0. How is this specifically implemented in the FreeFem++ code?

The issue is indeed to have an appropriate variational formulation, and to choose a good finite element space.
Here since you already set {\mathrm{div} B}=0 in \Omega, you don’t really need to impose {\mathrm{div} B}=0 in \partial\Omega because the latter follows from the former.
If you want nevertheless to enforce such condition you could add a term -\nabla\mathrm{div} B in the B equation. Then you have to take a space of the type P1 for B (you need \nabla B in L^2).
On the contrary if you don’t put such term you can take a space with only \nabla\times B\in L^2, such as RT03d or Edge03d, this should be more efficient.
I notice that you don’t have to really impose {\mathrm{div} B}=0 in \Omega because if it holds at t=0, the B equation implies that it will hold for all t>0 (but maybe not exactly at the discrete level).