Compute the slope of a boundary.

Question

Hello everybody. I would like to compute the slope of a boundary on a 2D FEniCS mesh.

For a simple example, consider a quadrilateral domain with vertices at (0,0), (2,1), (3,3), and (1,2). The bottom boundary is defined as the line between points (0,0) and (2,1). Now let the function $f(x)$ describe the vertical position of the bottom boundary as a function of $x$. In this case, $f(x) = \frac{1}{2}x$, and so ${\mathrm{d}f}/{\mathrm{d}x} = \frac{1}{2}$.

I need to compute ${\mathrm{d}f}/{\mathrm{d}x}$ along the boundary of a FEniCS mesh. The boundary won't necessarily be linear either - in fact, more than likely the boundary will be curved. The calculation needs to be done within FEniCS because some of the terms in the variational form depend on ${\mathrm{d}f}/{\mathrm{d}x}$.

Does anybody know a convenient way to perform a calculation like this? Thank you in advance.

Answer 1

Maybe you can use

n = FacetNormal(mesh)

Comments

I thought about using FaceNormal(mesh). But is there any way to print the nodal values of the normal vector though? And let's say I have a function g defined over the entire mesh:

# mesh and function space
mesh = UnitSquareMesh(3,3)
V = VectorFunctionSpace(mesh, "CG", 1)

# functions
g = Function(V)
n = FacetNormal(mesh)

Do you know of any way to populate g with values from the normal vector n?

---

no, the normal/slope is only defined on the boundary faces. Not on the nodes and not in the interior.

---

You can use this approach to get the a projection of the normal vector on vertices:

from dolfin import *
from mshr import *

domain = Sphere(Point(0, 0, 0), 1)
mesh = generate_mesh(domain, 50) 

V = VectorFunctionSpace(mesh, "CG", 1)

# Projection of the normal vector on P1 space
u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)

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

# Solve system
A = assemble(a, keep_diagonal=True)
b = assemble(L)
A.ident_zeros()

n = Function(V)
solve(A, n.vector(), b)
print assemble(div(n)*ds(domain=mesh))
plot(n, interactive=True)

---

That looks like a nice approach hernan_mella, I was looking for something like this.
One question though, it seems only normals on the exterior surface vertices are computed, all other nodes have vectors of length zero. How could this be adjusted?