Computing the boundary curvature of a 2D domain

206
views
1
6 months ago by
Hi all,

I am trying to obtain the curvature $\kappa$ of the boundary of a 2-dimensional domain by computing $\kappa=\nabla\cdot n$, where $n$ is the outer normal vector of the boundary. Since the normal $n$ obtained from the FacetNormal function is piecewise constant, I want to interpolate it to a function of a higher degree.

My current idea is to make Fenics solve the following problem:
$\int_\Gamma u\cdot v ds = \int_\Gamma n\cdot v, \quad \forall v \in H(div).$
I tested it on a unit circle, but did not get expected result. In particular, I expect $\int_\Gamma \nabla\cdot u ds=2\pi$, but the output of my code seems unbounded.

Any suggestions on the method? Or is there another way to get the curvature? Thank you.

from dolfin import *
from mshr import *

domain = Circle(Point(0, 0), 1)
mesh = generate_mesh(domain, 32)

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

u = TrialFunction(V)
v = TestFunction(V)
nn = FacetNormal(mesh)

a = inner(u, v)*ds
L = inner(nn, 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))
plot(n)

interactive()​
Community: FEniCS Project
I had the same problem the last year (see here). A workaround for this issue was using a measure defined on cells:

print(assemble(div(n)*dx))​

May be you can find a more elegant solution digging in the source code of FEMorph (take a look to the module SAD_geometry)

written 5 months ago by Hernán Mella