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