I am currently working on implementing the discontinuous Galerkin Finite Element Method (DG-FEM) for solving the Richards equation in a square domain with homogeneous Neumann boundary conditions on the lateral and bottom sides (q.n = 0). Initially, I assumed that these boundary conditions would naturally vanish in the weak form, but this doesn’t seem to be the case.

Below, I have outlined the weak form formulations for both the standard Finite Element Method (FEM) and the corresponding DG-FEM:

// RE: du/dt = div(Kh grad(h + z))

// Standard FEM:
int nn = 50;
mesh Th = square(nn, nn, [Lx, Ly]);
macro dn(u) (N.x * dx(u) + N.y * dy(u) ) // Define the normal derivative

I am encountering issues with the implementation and would greatly appreciate your assistance in resolving them. Any guidance or suggestions you can provide to correct this would be invaluable.

/______Continous Backward Euler _____________________/
macro grad(u) [dx(u),dy(u)] //
macro dn(u) (N’*grad(u) ) // def the normal derivative
real q0 = 0.5;
problem Richard(h, v, solver = GMRES) =
int2d(Th, qft = qf1pTlump)( Ah * h * v + Kh * (grad(h)'*grad(v)) )

int2d(Th, qft = qf1pTlump)( Ah * hm * v - (thetam - thetaold0) * v

Kh * dy(v) )

int1d(Th,3) (dt * q0 * v)
;
/_______________________________________________________________/
/________________Discontinous Galerkin BE _____________________/
problem Richard(h, v, solver = GMRES) =
int2d(Th, qft = qf1pTlump)( Ah * h * v + Kh * (grad(h)'*grad(v)) )

intalledges(Th)( ( jump(v) * mean(Kh * dn(h) )

pena * jump(h) * jump(v) ) / nTonEdge )

int2d(Th, qft = qf1pTlump)( Ah * hm * v - (thetam - thetaold0) * v

Kh * dy(v) )

int1d(Th,3) (dt * q0 * v)
;

I’m implementing the DG FEM corresponding to the Continuous FEM for the Richards equation. Please note that the Continuous FEM is working well. The flux in the Richards equation is
q = -Kh∇(h+y), where h is the unknown.

In the DG FEM implementation, I need to ensure that the no-flux (homogeneous Neumann) boundary conditions are correctly applied. Specifically, when setting intalledges, I want to ensure it considers all edges except those with no-flux boundary conditions.

Can you provide guidance on how to exclude contributions from the edges with no-flux boundary conditions while implementing the DG FEM in FreeFem++?