Compute flux: how to get the unit normal vector?

Question

Hi,
I am a new user of Fenics, and i don't understand how to compute a flux...
I would like to calculate integral(grad u.n ds), knowing grad u.
My problem is I don't get it how to define n (the unit normal vector)

Here is a simple version of the code I use

# Import libraries
from dolfin import *
# Create the mesh 
mesh = UnitCubeMesh(4,4,4)
# Discrete function space V
V = FunctionSpace(mesh, 'Lagrange', 1)
# Boundary condition
tolerance_boundary_location = 1e-15
def apply_boundary_omega1(x):
    return ( abs(x[0])  < tolerance_boundary_location or abs(x[0]-1)  < tolerance_boundary_location
            or abs(x[1])  < tolerance_boundary_location or abs(x[1]-1)  < tolerance_boundary_location
            or abs(x[2])  < tolerance_boundary_location or abs(x[2]-1)  < tolerance_boundary_location )
bc_value = Expression('-400*x[0] + 500')
bc = DirichletBC(V, bc_value, apply_boundary_omega1)
# Define trial function u
u_trial = TrialFunction(V)
# Define test functions v
v_test = TestFunction(V)
# Define function f
f = Constant(0)
# Integrands of a(u,v)
a = inner (nabla_grad(u_trial),nabla_grad(v_test))*dx
# Integrands of L(v)
L = f*v_test*dx
# Define the unknown function to be computed
u = Function(V)
# Solve the PDE variational problem
solve(a == L, u, bc)
# Plot the solution
plot(u, interactive=True)
# Calculate the gradient through the projection method, using built-in function project for the calculation
gradient_u = project(-grad(u),VectorFunctionSpace(mesh, 'Lagrange', 1))
# Plot the gradient
plot(gradient_u, interactive=True)

Then, what is the best solution to compute integral (grad u .n ds) on the the plane x=1 and normal to it (then n=[0, 1])
- should i use the function FacetNormal ? (this function is not clear for me)
- should i define manually n such as n=[0, 1] and assemble ? (it's not working)
- should i perform a loop on each cell, grab each local grad u, calculate locally grad u.n and sum all the values ? (seems not a very elegant solution...)

Any help will be very appreciate.
Regards,

Answer 1

Hi,
definitely the two first methods are the way to do it. A simple manner to do it is the next:

from dolfin import *

mesh = UnitCubeMesh(4,4,4)
F = FunctionSpace(mesh, "CG", 2)
V = VectorFunctionSpace(mesh, "CG", 2, dim=3)

u = interpolate(Expression("x[0]*x[1]*x[2]*sin(x[0]*x[1]*x[2])"), F)
grad_u = project(grad(u), V)

# Create subdomain (x0 = 1)
class Plane(SubDomain):
  def inside(self, x, on_boundary):
    return x[0] > 1.0 - DOLFIN_EPS

# Mark facets
facets = FacetFunction("size_t", mesh)
Plane().mark(facets, 1)
ds = Measure("ds")[facets]

### First method ###
# Define facet normal vector (built-in method)
n = FacetNormal(mesh)
flux_1 = assemble(dot(grad_u, n)*ds(1))

### Second method ###
# Manually define the normal vector
n = Constant((1.0,0.0,0.0))
flux_2 = assemble(dot(grad_u, n)*ds(1))

print "flux 1: ", flux_1
print "flux 2: ", flux_2

Here you can find the meaning of FacetNormal for a mesh of triangles in 2D or 3D.

Comments

Thanks you, it works perfectly.
If i may ask, there are 2 points I would like to discuss more:

The line: ds = Measure("ds")[facets] is working but generates a warning:
Notation dx[meshfunction] is deprecated. Please use dx(subdomain_data=meshfunction) instead.
Then, what is the new syntax ?
I think it is ds = Measure("ds")(*=meshfunction) but i don't know what should be *

I can plot the plane with the following line to check my subdomain
plot(facets, interactive=True)
Then, what would be great is to plot the solution u only on this subdomain
My guess is I have to define a new function space
V_plane = FunctionSpace(facets, "CG", 2) # It doesnot work because facets is not a mesh
And then, 'project' the solution on this new function space

Regards,

---

The updated syntax would be

ds = ds(subdomain_data = meshfunction)

See here for more information about your second question (maybe you will need to update the code to the actual version of dolfin).