I am trying to add a surface tension-like term to my equations as a Von-Neumann boundary condition. This requires the curvature of my (2D) surface to be calculated. I’ve seen that this has been done in Dolfin (How to compute curvature of a boundary) but is there an equivalent for Dolfinx?

As an attempt, I’ve used:

Facet normal and facet tangent vector approximation in dolfinx · GitHub

Normal in dolfinx

to project the normal into a vector space (Discontinuous) and then tried to use this function in my equation but this didn’t seem to work.

Alternatively, I might have misunderstood this but from this post Extract normal vector data at boundary in dolfinx - #7 by dokken, using w*ufl.div(ufl.FacetNormal)*ufl.ds (with w a testfunction) should now be possible?

Thanks

If you are using straight edges, the divergence of a facetnormal is zero. See for instance: Normal vector in manifold - #6 by sdcardozof

Thank you for this! I had one quick follow up question which is that as my mesh is changing, I want to ensure the relevant values going into the weak formulation update as the mesh changes. Does this mean I need to use a fem.Function or can I put “n_oriented” directly into my weak formulation and it will do this?

If the mesh changes you would have to interpolate into `n_oriented`

again.

To ensure that I understand your question it would ne beneficial with a minimal reproducible example.

So as background, I am solving a Stokes equation (with a mixed function space formulation) with elliptical/circular mesh denoted by “domain”.

```
n = ufl.FacetNormal(domain)
x = ufl.SpatialCoordinate(domain)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n
# Constant in surface curvature term
gamma_i = fem.Constant(domain, PETSc.ScalarType(gamma))
# Term of interest which is a Neumann BC with it involving surface tension
# and added to "RHS" of Stokes equation
L_v += fem.form(ufl.inner(gamma_i*ufl.div(n_oriented)*n,v)*ufl.ds)
# Omitted specifying solver, creating matrix as well as defining function space to shift mesh nodes/points
# Defining vector which will contain result
U_out = fem.Function(V_stokes)
# Loop
for i in range(0,100):
# Removed other parts in solving (assembling vector, matrix etc.)
ksp_v.solve(bvec_v, U_out.vector)
u_out = U_out.sub(0).collapse()
# Shift mesh points by u_final (Extra code to make shifting the domain work correctly removed)
domain.geometry.x[:,:domain.geometry.dim] = u_final.x.array.reshape((-1,
domain.geometry.dim)) + domain.geometry.x[:,:domain.geometry.dim]
```

u_final is u_out * dt [Time step] in the “domain” function space which is displacement for mesh points.

I am using the loop because as the mesh is displaced I need to remesh/ restore its “quality” after a certain amount of time steps. So from my understanding, each time I call solve, it finds a and L again for the new adjusted domain/mesh with the nodes at different positions. Does this happen for “n_oriented” too or is it finding “n_oriented” once at the start and using the same value throughout the entire loop?

Thanks and I am happy to add more if my code is too sparse

As `n_oriented`

is a pure ufl expression, it gets updated with any change in the mesh geometry (as long as you call assemble inside your loop).