### How to compute boundary flux?

212
views
1
3 months ago by
Dear all,

I'd like to compute the boundary flux of a function: $\frac{\partial u}{\partial\vec n}$, where $\vec n$ is the normal vector at the boundary.
I need values for every DOF, and NOT the total flux, i.e., the integral over all the boundary,

Ideally, something like the following would work:
# V = FunctionSpace(...)
n = FacetNormal(mesh)
project(dot(grad(u), n), V)​

However, it does not and I haven't found a simple way of doing what seems a very standard operation to me?

If not I could iterate over the facets and compute dot(grad(u), n) by hand, or solve $\int_{\partial\Omega} wz = \int_{\partial\Omega} \nabla u\cdot\vec n z$ via least squares (or eliminating all interior DOFs with DirichletBCs). But this seems overly complicated.

Thanks!

Community: FEniCS Project

0
3 months ago by
The code you posted will not work, and should not work, because the outward-pointing-unit-normal vector is not defined in the interior.  That said, a method that worked for me in the past is to project the facet normals defined by the geometry to a function, then solve the dot product with grad(u) using standard projection

https://fenicsproject.org/qa/8931/get-normal-component-of-a-3d-vector-field

Try that

edit: to clarify from that old post, here:

from fenics import *

mesh    = UnitCubeMesh(20,20,20)
V       = VectorFunctionSpace(mesh, 'CG', 1)
n       = FacetNormal(mesh)

u       = TrialFunction(V)
v       = TestFunction(V)

a = inner(u,v) * ds
l = inner(n,v) * ds

A = assemble(a, keep_diagonal=True)
A.ident_zeros()
b = assemble(l)

n_f = Function(V)

solve(A, n_f.vector(), b)

File('n.pvd') << n_f

note that you're going to get some noisy crap around the edges because the normal vector at vertices will be averages of the normals defined at adjacent cell faces; this is the reason why integrated values of flux are more accurate and useful.  If you have a surface that is defined implicitly, you have an exact equation for the surface, and can use this in place of the FEniCS-computed "n".  The normal you get from projection in either case should converge to the correct normal with mesh refinement.
1
A.ident_zeros() is nicer than using DirichletBCs as I wrote in my question. Thanks for the hint
written 3 months ago by David