Adding Surface Curvature Terms

dolfinx
1

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

2

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

1 Like
3

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?

4

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.

5

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

6

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).

1 Like