Hello all. I have a question with regards to how FEniCS is calculating the gradient of a scalar function at the boundary. I am trying to solve the Poisson equation where my source term f is a vector function with discrete values implemented on piecewise constant DG function space as such:
f = Function(V_pwc)
f_values = f.vector().get_local()
f_values[cellno] = 5e7
f.vector().set_local(f_values)
The u function space is 1st order lagrangian. On a 2D unit-square grid, the dirichlet BCs are all set to 1. If I set an internal element to the source value (see picture below), I can calculate the ratio of the flow in due to the source, and the flow out of the boundaries, expectign a ratio of approximately 1, which I do get.
flux_due_to_source = element_volume*f_value
flux_through_boundaries = dot(-grad(u), n)*ds
However, if the source is at a boundary (such as in the image below), the ratio of the fluxes is now 3
And if the source only has a node on the boundary, not an edge, the ratio of the fluxes is 1.5.
I have been banging my head against a wall trying to figure out what I’m doing wrong. Is there fundamentally a problem with defining an elemental source at a dirichlet boundary?