Computing the boundary curvature of a 2D domain
6 months ago by
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.