### Computing the boundary curvature of a 2D domain

206

views

1

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.

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

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

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