Take derivative along 1D mesh

dolfinx
1

Hi,

I am working on a 1D problem embedded in 3D space.

The problem requires taking the derivative of a function in the direction of the mesh. I have an inelgant solution by projecting the 3D derivative (although I’m unsure how this is defined as mesh.topology.dim == 1) along the direction of the mesh, but this seems to be quite inelegant and potentially inaccurate.

Is there a better way of doing this? I understand that the integral equivalent of this problem can be solved neatly using ufl.Measure(). Is there a similarly neat solution for differentiation?

MWE:

# %% Imports
import gmsh
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
from dolfinx import io
from dolfinx.fem.petsc import LinearProblem
import ufl
import ufl.constant
import numpy as np

# Seems to work but is not an elegant solution for a mesh which is topologically 1D

def build_mesh():
    """create a 1D mesh that lives in 3D space"""
    gmsh.initialize()
    gmsh.option.setNumber('General.Verbosity', 1)

    pts = []
    lines = []
    meshSize = 0.1
    pts.append(gmsh.model.geo.addPoint(0., 0., 0., meshSize=meshSize, tag=1))
    pts.append(gmsh.model.geo.addPoint(1., 1., 0., meshSize=meshSize, tag=2))
    lines.append(gmsh.model.geo.addLine(1, 2))
    gmsh.model.addPhysicalGroup(dim=1, tags = lines)
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(1)
    (msh, _, _) = io.gmshio.model_to_mesh(gmsh.model, comm=MPI.COMM_WORLD, rank=0)

    gmsh.finalize()
    return msh # returns a mesh with cells along the vector [1, 1, 0]


msh = build_mesh() # create 1D mesh
print("mesh topological dimension = ", msh.topology.dim, " mesh geometric dimension = ", msh.geometry.dim)

V = fem.functionspace(msh, ("DG", 1)) # create 1D functionspace
U = fem.Function(V)

U.interpolate(lambda x: x[0] + x[1]) # set initial condition
l = fem.Function(V)
l.interpolate(lambda x: x[0] + x[1] + x[2]) # set mesh vertex position function
# ufl.grad(l) seems to return l.dx(0) + l.dx(1) + l.dx(2) despite the 1d topological dim
l_hat = ufl.grad(l)/ ufl.sqrt(ufl.dot(ufl.grad(l), ufl.grad(l))) #create normalised mesh direction vector

u_trial = ufl.TrialFunction(V)
v_test = ufl.TestFunction(V)
dx = ufl.Measure("dx", domain=msh)

a_expr = ufl.inner(u_trial, v_test)*dx
L_expr = ufl.inner(ufl.dot(ufl.grad(U), l_hat), v_test)*dx #the variational form of the projection of grad(U) in the mesh direction. e.g. dUdl

problem = LinearProblem(a_expr, L_expr)
dUdl = problem.solve()
dUdl = dUdl.x.array

assert np.allclose(dUdl, np.sqrt(2)), "test_failed"
print("test passed")
2

You can take partial derivatives if you prefer:

L_expr = ufl.inner(np.sqrt(2)/2*(ufl.Dx(U, 0) + ufl.Dx(U, 1)), v_test)*dx

EDIT
I found this solution after a little bit of digging, it’s more elegant:

verts = ufl.classes.CellVertices(msh)
tangent = ufl.as_vector([verts[1, i] - verts[0, i] for i in range(3)])
utangent = tangent / ufl.sqrt(ufl.dot(tangent, tangent))
L_expr = ufl.inner(ufl.dot(ufl.grad(U), utangent), v_test)*dx #the variational form of the projection of grad(U) in the mesh direction. e.g. dUdl