Interpolate() function not working properly with BDM elements

1

Hello, I am using FEniCS 2019.01 and am experiencing and issue when trying to interpolate a specific Expression to the BDM1 space. The problem is set on a isosceles triangle element with vertices (1,0), (0,1), (-1,0). The considered function is v=(0,x^2). As you can see from the output of the code below, the interpolated function seems to be Iv = (\frac{2x}{3}, -\frac{y}{3}+\frac{1}{9}).

Yet when calculating the interpolated function by hand and in NGSolve with the method provided in https://ngsolve.org/docu/latest/i-tutorials/unit-2.10-dualbasis/dualbasis.html, I get to the interpolant Iv = (\frac{x}{2},-\frac{y}{2}+\frac{1}{3}).

So is there an issue here in FEniCS, or am I making a mistake?

from dolfin import *

editor = MeshEditor()
mesh = Mesh()
editor.open(mesh, 'triangle', 2, 2)
editor.init_vertices(3)
editor.add_vertex(0,[-1,0])
editor.add_vertex(1,[1,0])
editor.add_vertex(2,[0,1])
editor.init_cells(1)
editor.add_cell(0, [0,1,2])
editor.close()
mesh.init()
x = SpatialCoordinate(mesh)

W_BDM = FunctionSpace(mesh, "BDM", 1)

v = Expression(('0','x[0]*x[0]'), degree = 2)
Iv_guess = as_vector([2.0*x[0]/3.0,-x[1]/3.0+1.0/9.0])

Iv = interpolate(v, W_BDM)

print('L2-norm of Iv-Iv_guess:', str(sqrt(assemble((Iv-Iv_guess)**2 * dx(mesh)))))

output:

L2-norm of Iv-Iv_guess: 2.127340761518145e-16

2

After more reflection, I really think there is something wrong with the BDM interpolation using the interpolate() function:

The BDM interpolation operator is supposed to map divergence-free functions to divergence-free BDM functions. As you can see in the above example, the function which FEniCS calculates is clearly not divergence-free even though the input function is.

Is my observation correct and is there a way to fix this?

3

I did some more experimentation, if anybody is interested in a follow up:

Using the ideas from
https://fenicsproject.org/qa/5245/nedelec-interpolation/
I got to the following program to “manually” calculate the BDM-interpolant of a function on a triangle, which seems to be working very nice, as it reproduces the results I get from the NGSolve interpolation.

Anyhow, this still leaves the question what the internal interpolate() function does, and why it seems to do the wrong thing when interpolating to a BDM function space.

from dolfin import *
import numpy as np

order = 2

editor = MeshEditor()
mesh = Mesh()
editor.open(mesh, 'triangle', 2, 2)
editor.init_vertices(4)
editor.add_vertex(0,[0,0])
editor.add_vertex(1,[1,0])
editor.add_vertex(2,[0,1])
editor.init_cells(1)
editor.add_cell(0, [0,1,2])
editor.close()
mesh.init()

def BDMi_Ned_dofs(u0, mesh, r):
    V = FunctionSpace(mesh, 'BDM', r)
    u = TrialFunction(V)
    v = TestFunction(V)
    if r == 1:
        W = FunctionSpace(mesh, 'BDM', 1)
        v0 = TestFunction(W)
    else:
        V0 = FiniteElement('BDM', mesh.ufl_cell(), r)
        V1 = FiniteElement('N1curl', mesh.ufl_cell(), r-1)
        W = FunctionSpace(mesh, MixedElement([V0,V1])) 
        v0, v1 = TestFunctions(W)
        
    n = FacetNormal(mesh)
    a = (inner(u, n) * inner(v0, n))('+') * dS  +  inner(u, n) * inner(v0, n) * ds
    L = (inner(u0, n) * inner(v0, n))('+') * dS  +  inner(u0, n) * inner(v0, n) * ds
    if r > 1:
        a = a + inner(u, v1) * dx
        L = L + inner(u0, v1) * dx
        
    A = assemble(a)
    A_reduced = assemble(inner(u,v)*dx(mesh))
    b = assemble(L)
    b_reduced = assemble(inner(v,Function(V))*dx(mesh))
    
    n_facet_dofs = A.size(0) - 2*(A.size(0)-A.size(1))
    n_internal_dofs = V.dim() - n_facet_dofs
    cols = np.array([i for i in range(A.size(1))], dtype=np.intc)
    
    block = np.ndarray((V.dim(),V.dim()), dtype=np.float_)
    for i in range(n_facet_dofs):
        block[i,:] = A.getrow(i)[1][:]
        b_reduced[i] = b[i]
    for i in range(n_internal_dofs):
        block[n_facet_dofs+i,:] = A.getrow(n_facet_dofs+n_internal_dofs+i)[1][:]
        b_reduced[n_facet_dofs+i] = b[n_facet_dofs+n_internal_dofs+i]
    
    A_reduced.set(block, cols, cols)
    A_reduced.apply('insert')
    b_reduced.apply('insert')
    
    u = Function(V)
    solve(A_reduced, u.vector(), b_reduced)
    return u

V = FunctionSpace(mesh, 'BDM', order)
g = Expression(('0','pow(x[0],3)'), degree = 3)

#u = interpolate(g,V)
u = BDMi_Ned_dofs(g, mesh, order)
print('Interpolation manual, Nedelec dofs:')
print('\t L2 norm u[0]', assemble(u[0]*u[0]*dx(mesh)))
print('\t L2 norm u[1]', assemble(u[1]*u[1]*dx(mesh)))
print('\t L2 norm div(u)', assemble(div(u)*div(u)*dx(mesh)))
print('\n')

u = interpolate(g, V)
print('Interpolation automatic:')
print('\t L2 norm u[0]', assemble(u[0]*u[0]*dx(mesh)))
print('\t L2 norm u[1]', assemble(u[1]*u[1]*dx(mesh)))
print('\t L2 norm div(u)', assemble(div(u)*div(u)*dx(mesh)))
1 Like
4

:orange_heart: :orange_heart:

5

Is this problem still present in FenicsX?

6

No, it shouldn’t be present any longer, as we have completely rewritten how interpolation is done.