Curvature by finite elements

variational formulation
1

Hello everyone.

I want to compute the curvature of a closed curve by finite elements.

The variational problem reads:

image
image925×150 20.7 KB

So I tried to formulate it in FEniCS as:

from fenics import *

mesh = Mesh('CircleSurface.xml')

V = VectorFunctionSpace(mesh, 'P', 2)

H = TrialFunction(V)
phi = TestFunction(V)

def Identity(u):
    return Expression(("x[0]","x[1]"), degree = 1)

Idt = project(Identity(H), V)

F = inner(phi, H)*dx - inner(nabla_grad(phi), nabla_grad(Idt))*dx
a, L = lhs(F), rhs(F)
H = Function(V)
solve(a==L, H)

However, the solution is NaN at every point. I checked and Idt is also NaN. I really do not know why. As far as I know, the expression I used is correct, why is it NaN?

In addition, I would like to know if there is any source of tangential gradient, so far I am using the flat gradient assuming that the integration takes into account the surface in the Jacobian.

I implemented this problem in matlab and I used a map from the curved element to a reference flat 1D domain, there the Jacobian takes into account the fact that the domain is a surface.

Is there any similar approach in FEniCS?

I also want to add a short view of the mesh: (‘CircleSurface.xml’)

<?xml version="1.0"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <mesh celltype="interval" dim="2">
    <vertices size="100">
      <vertex index="0" x="1.0000000000000000e+00" y="0.0000000000000000e+00"/>
      ...
      <vertex index="99" x="9.9802672842827156e-01" y="-6.2790519529313263e-02"/>
    </vertices>
    <cells size="100">
      <interval index="0" v0="0" v1="0" />
      ...
      <interval index="99" v0="99" v1="99" />
    </cells>
    <data />
  </mesh>
</dolfin>
2

I could solve the problem. It seems that project does not work well as indicated here so I used interpolate instead. In addition I defined the identity as an expression-although I am not sure if that was also affecting the code-as:

Identity = Expression(("x[0]","x[1]"), degree=1)
Idt = interpolate(Identity, V)

Finally, here Rogens indicates that the approximation only works with linear elements, so I changed the degree of the vector space to 1 and it worked. I am going to run some extra simulations to see how it goes with more complex geometries.

I also want to add that the fact that the implementation only works with linear elements is a big issue in this particular case. Heine showed that the order of convergence of this isoparametric approach is degree-1 (degree of Lagrange elements). So in the case of linear elements, it is not convergent.

If anyone finds a way to use higher-order elements I would really appreciate if they share it with me. I think this issue can be overcome if we could import a quadratic mesh including the mid nodes of each element. I guess it might be already possible but I have not found how to do it.

3

Hi, @davhernandezari.

Have you found an answer to this?

I ask that since I’m going to work with diffusion in manifolds, that way I was wondering if is necessary to use the correspondent of the surface gradient \nabla_S = (\mathbf{I} - \mathbf{n}\otimes \mathbf{n})\nabla, or that grad is equivalent to \nabla_S in this situation.

4

Hi @fforteneto.

From the work by Rogens it is equivalent but only with linear elements. I have not tried it though.

1 Like