Interpolate into complex-valued vector function

I was trying to test a complex-valued PDE on a standard rectangular domain and noticed some issues with the way I was interpolating and expression into a function.

Here is a MWE to get a function \boldsymbol{u} = \exp(-\mathrm{i} \boldsymbol{k} \cdot \boldsymbol{x}) \boldsymbol{k} for a given plane-wave direction vector, \boldsymbol{k}.

import numpy as np
from mpi4py import MPI
from dolfinx.mesh import create_rectangle
from dolfinx.fem import FunctionSpace, Function

msh = create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 1.0)), (20, 20))
msh.name = 'dolfinx_mesh'

def foo(x):
    k = [1.0, 0.0]
    return np.vstack((
        np.exp(-1j*(k[0]*x[0]+k[1]*x[1])) * k[0],
        np.exp(-1j*(k[0]*x[0]+k[1]*x[1])) * k[1]
    ))

V = FunctionSpace(mesh, ('Raviart-Thomas', 1))

u = Function(V, dtype=np.complex128)
u.interpolate(foo)

Results in some issues along the Y-coordinate (which is supposed to be 0. since it is multiplied explicitly with it)

This bar on the left-hand side that you notice is observed also at other boundaries (with more complex meshes). The width of this bar does reduce with increased RT-order and mesh refinement, but not entirely and artifacts do remain in a few cells.

Could someone point me to something I am doing wrong in interpolating a vector-valued expression to a Raviart-Thomas function space? I borrowed this interpolation template from test_interpolation.py.

1 Like

I don’t think you’re doing anything wrong. You should, however, consider interpolating to a DG-1 function space and using VTXWriter for exporting results, as noted here:

Below, I compare the real part of the y-component, exported with interpolation into a DG-1 space (plots labeled “VTK”) and without interpolation into a DG-1 space (plots labeled “XDMF”). You can see that the error depends on the type of discretization (left, right, or crossed).

In general, the erroneous nonzero y-component is very small, as you can see by plotting your interpolated field with a “Glyph” filter, so I would guess that the interpolation is giving you the best representation of the plane wave that you can achieve with RT elements.

3 Likes

I see, thank you! I could reproduce the plots that you got once I projected the function into DG space.

Just as a remark, this discretization error (although a couple orders small) still seems significant since they appear when complex valued functions are utilized. I had not encountered these earlier when I used two-part real RT space to represent the complex functions (in dolfin, although I believe it did some post-processing while writing RT into XDMF).

Depending on what function you used with XDMFFile,

  • write(u) (interpolates into a CG-1 space)
  • write_checkpoint(u) (for some orders represented with: XDMF Model and Format - XdmfWeb section finite element function).