Indicator function to implement boundary conditions

General
1

Hello, attached I have an implementation of an HDG method for a Poisson problem.
bc_problem_poisson.py (2.2 KB)

Enforcing Neumann boundary conditions with the below code gives the expected convergence rates:

ds = dx(element_boundary = True)
f += ((grad_uex * n) * vbar.Trace()) * ds(definedon=mesh.Boundaries('Neumann'))

Now instead of integrating over ds(definedon=mesh.Boundaries(‘Neumann’)), I would like to implement the Neumann boundary conditions by integrating over all element boundaries, and using an indicator function (or similar) so that the contribution from internal element boundaries is 0.
My attempt:

chiN = GridFunction(FacetFESpace(mesh, order=0, dirichlet='Neumann'))
chiN.Set(CF(1.0), definedon=mesh.Boundaries('Neumann'))
f += ((grad_uex * n) * vbar * chiN) * ds

This does not converge. How can I modify this to achieve the same result as I get when integrating over ds(definedon=mesh.Boundaries(‘Neumann’))?

Thanks,
Elizabeth

2

you have to reset the chiN after refinement:

for i in range(4):
    mesh.Refine()
    chiN.Set(CF(1.0), definedon=mesh.Boundaries('Neumann'))
    Solve()

then I got the output

 DOFs     E      r
 25808 4.70e-07 4.00
102816 2.93e-08 4.00
410432 1.83e-09 4.00

It is relatively new that spaces and gridfunctions are updated after mesh refinement, but FacetFESpace has no automatic prolongation (since it is not clear what that should be).

3

That works for me, thanks!