### How to compute boundary flux?

502

views

1

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:

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

### 1 Answer

0

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

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
9 months ago by
David Nolte

Please login to add an answer/comment or follow this question.