DOLFIN-X: Curl elements treated differently from DOLFIN

1

Hi,

I am trying to solve the vector potential formulation of Maxwell’s equations in DOLFIN-X, which is preferred for its complex number support. To do so, I’m using 1st-order HCurl elements. However, it seems that DOLFIN-X handles these elements differently from how DOLFIN does, which is problematic because of the lack of documentation for DOLFIN-X. As an example, I’ve produced the following code using DOLFIN:

from dolfin import *

#Defines mesh and function space
mesh = UnitCubeMesh(1,1,1)
V = FunctionSpace(mesh, 'N1curl', 1)

#Defines boundary
def boundary(x, on_boundary):
    return on_boundary

#Defines boundary function
u_expr = Expression(('1.0','0.0','0.0'), element = V.ufl_element())
u_func = interpolate(u_expr, V)
bc = DirichletBC(V, u_func, boundary)

#Writes out file
outfile = XDMFFile(mesh.mpi_comm(), 'BC_Dolfin.xdmf')
outfile.write(u_func)
outfile.close()

#Obtains test and trial funtions
u = TestFunction(V)
v = TrialFunction(V)

#Defines variational problem
f = Function(V)
a = inner(curl(v), curl(u)) * dx
L = -inner(v, f) * dx

#Solves variation problem
u = Function(V)
solve(a == L, u, bc)

#Writes out results
outfile = XDMFFile(mesh.mpi_comm(), 'Results_Dolfin.xdmf')
outfile.write(u)
outfile.close()

An equivalent script in DOLFIN-X (I think) is as follows:

from dolfin import *
from dolfin.io import XDMFFile
from ufl import curl, inner, dx
import numpy as np

#Defines mesh and function space
mesh = UnitCubeMesh(MPI.comm_world, 1, 1, 1)
V = FunctionSpace(mesh, ('N1curl', 1))

#Defines boundary
def boundary(x, on_boundary):
    #Hacky workaround: fails if only on_boundary is used
    return np.logical_or(on_boundary, x[:, 0] == 1e3)

#Defines boundary function
@function.expression.numba_eval
def u_bc(values, x, cell):
    values[:, 0] = 1.0
    values[:, 1] = 0.0
    values[:, 2] = 0.0
u_expr = Expression(u_bc, shape=(3,))
u_func = interpolate(u_expr, V)
bc = DirichletBC(V, u_func, boundary)

#Writes out file
outfile = XDMFFile(MPI.comm_world, 'BC_DolfinX.xdmf')
outfile.write(u_func)
outfile.close()

#Obtains test and trial funtions
u = TestFunction(V)
v = TrialFunction(V)

#Defines variational problem
f = Function(V)
a = inner(curl(v), curl(u)) * dx
L = -inner(v, f) * dx

#Solves variation problem
u = Function(V)
solve(a == L, u, bc)

#Writes out results
outfile = XDMFFile(MPI.comm_world, 'Results_DolfinX.xdmf')
outfile.write(u)
outfile.close()

These scripts produce two XDMF files for visualization in paraview: one for the boundary value function and the other for the solution. Shown below, these two scripts produce very different boundary functions.

The left image is the output from DOLFIN and the right is the output from DOLFIN-X

03%20AM
Screen Shot 2019-06-12 at 11.48.03 AM.jpg2220×942 139 KB

The results become even stranger if a spherical mesh is used. However, they produce the same results and the same (presumably incorrect) solution if Lagrange (CG) vector elements are used instead.

Is there some key difference in how DOLFIN-X handles N1curl elements compared to DOLFIN? Do the degrees of freedom in the boundary function need to be treated differently? I have limited experience with FEM and am unsure of how to properly work with these elements to solve my problem, so any help would be appreciated. Thanks!

2 Likes