# Inforce homogeneous Neumann BCs using DG FEM

Dear experts,

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

varf Richard(h, v, solver = UMFPACK) =
int2d(Th, qft = qf1pTlump)( Ah * h * v + Kh * (dx(h) * dx(v) + dy(h) * dy(v)) )

• int2d(Th, qft = qf1pTlump)( Ah * hm * v - (thetam - thetaold0) * v
• Kh * dy(v) )
• int1d(Th, 3)(q0 * v)
;

// DG-FEM:
problem Richard(h, v, solver = UMFPACK) =
int2d(Th, qft = qf1pTlump)( Ah * h * v + Kh * (dy(h) * dy(v) + dx(h) * dx(v)) )
+ intalledges(Th)( ( jump(v) * mean(Kh * dn(h)) - jump(h) * mean(Kh * dn(v))
+ pena * jump(h) * jump(v) ) / nTonEdge )
- int2d(Th, qft = qf1pTlump)( Ah * hm * v - (thetanew - thetaold0) * v - Kh * dy(v) )
- int1d(Th, 2)(0.15 * dn(v) + pena * 0.15 * v)
;

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.

Best regards,

Noureddine

Sorry, DG FEM is not well defined, This a a lot of formulation , so what the formulation.

/______Continous Backward Euler _____________________/
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++?

You can just multiply by `(nTonEdge-1)` to delete the boundary edges,
see About intalledges and internal edges - #2 by fb77