Shape optimization question: div(n)?

Question

Dear all,

Makes sense to use in a weak form (defined only on boundary) the term div(n) (with n = FacetNormal(mesh))?. Or it's better project the normal vector to a function space defined over the facets? (something like Raviart-Thomas)

Thanks in advance!

Comments

Hi, FacetNormal is a facet-wise constant so its div is zero. As a consequence a sphere has zero mean curvature in the following computation

from dolfin import *
from mshr import *

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

n = FacetNormal(mesh)
print assemble(div(n)*ds(domain=mesh))

---

Thanks for your answer. I was expecting that FEniCS would allow to me "magically" perform this operation. Anyway, I was able to solve my problem with this:

from dolfin import *
from mshr import *

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

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

# Projection of the normal vector on P2 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)

---

Dear Miroslav,

Do you know why using my approach I'm not getting

$$\int_{\partial\Omega} \text{div}~n~ds = 4\pi?$$

On the other hand, using assemble(div(n)*dx(domain=mesh)) the result is $4\pi$ (the result that I was expecting with the previous integral operation).

---

Hi, not at the moment but I will have a look. Interesting, the quantity assemble(div(n)*ds(domain=mesh)) is unbounded.