Inforce homogeneous Neumann BCs using DG FEM

General Discussion
1

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

2

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

3

/______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++?

4

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