Interpolate from enriched space to non-enriched space with identical meshes.

Question

I know that there are some difficulties interpolating from Enriched function spaces to non-enriched spaces.

The attached script here used to work, but now the "stabilized" solution is not correct. Is there a way to fix this with dolfin 1.6?

Thanks again,
Evan

from fenics import *

mesh = IntervalMesh(10,0,1)

Q    = FunctionSpace(mesh, 'CG', 1)
B    = FunctionSpace(mesh, 'B',  2)
M    = Q + B

def left(x, on_boundary):
  return on_boundary and x[0] == 0

leftBC  = DirichletBC(Q, 1.0, left)

kappa = 1.0/500.0
s     = 1.0
f     = Constant(0.0)

#==============================================================================
# unstabilized :

u   = TrialFunction(Q)
v   = TestFunction(Q)
us  = Function(Q)
uf1 = Function(Q)

a  = + kappa * u.dx(0) * v.dx(0) * dx \
     + s * u * v * dx
L  = f * v * dx

solve(a == L, us, leftBC)

uf1.interpolate(us)

#==============================================================================
# stabilized :

u   = TrialFunction(M)
v   = TestFunction(M)
us  = Function(M)
uf2 = Function(Q)

a  = + kappa * u.dx(0) * v.dx(0) * dx \
     + s * u * v * dx
L  = f * v * dx

solve(a == L, us, leftBC)

uf2.interpolate(us)

#==============================================================================
# plotting :
from pylab import *

t    = mesh.coordinates()[:,0][::-1]
uf1  = uf1.vector().array()
uf2  = uf2.vector().array()

x    = linspace(0, 1, 1000)
ue   = (exp(-10*sqrt(5)*(x-2)) + exp(10*sqrt(5)*x))/(1 + exp(20*sqrt(5)))

plot(t, uf1, 'r',   lw=2.0, label=r"$u_1$")
plot(t, uf2, 'k',   lw=2.0, label=r"$u_2$")
plot(x, ue,  'k--', lw=2.0, label=r"$u_a$")

xlabel(r'$x$')
ylabel(r'$u$')
legend()
xlim([0,0.3])
grid()
show()

Answer 1

Interpolation from enriched to non-enriched on the same mesh should be fine. Watch https://bitbucket.org/fenics-project/ffc/issues/69 for the resolution of problems with enriched element. Any known problem with enriched element should raise an error rather than compute a rubbish. If not it should be reported.

Anyway your code contains a mistake. You should rather do:

uD = project(Constant(1.0), M) # workaround FFC issue #69
leftBC  = DirichletBC(M, uD, left)

Comments

First, I upgraded to v. 2016.1.0 recently. I tried your workaround, but am getting the same rubbish output as before.

I have the full code with a plot to demonstrate:

from fenics import *

mesh = IntervalMesh(10,0,1)

Q    = FunctionSpace(mesh, 'CG', 1)
B    = FunctionSpace(mesh, 'B',  2)
M    = Q + B

def left(x, on_boundary):
  return on_boundary and x[0] == 0

leftBC  = DirichletBC(Q, 1.0, left)

kappa = Constant(1.0/500.0)
s     = Constant(1.0)
f     = Constant(0.0)

#==============================================================================
# standard Galerkin solution :

u   = TrialFunction(Q)
v   = TestFunction(Q)
us  = Function(Q)
uf1 = Function(Q)

a  = + kappa * u.dx(0) * v.dx(0) * dx \
     + s * u * v * dx
L  = f * v * dx

solve(a == L, us, leftBC)

uf1.interpolate(us)

#==============================================================================
# enriched space solution :

uD      = project(Constant(1.0), M)
leftBC  = DirichletBC(Q, uD, left)

u   = TrialFunction(M)
v   = TestFunction(M)
us  = Function(M)
uf2 = Function(Q)

a  = + kappa * u.dx(0) * v.dx(0) * dx \
     + s * u * v * dx
L  = f * v * dx

solve(a == L, us, leftBC)

uf2.interpolate(us)
#==============================================================================  
# plotting :
from pylab import *

t    = mesh.coordinates()[:,0][::-1]
uf1  = uf1.vector().array()
uf2  = uf2.vector().array()

x    = linspace(0, 1, 1000)
ue   = (exp(-10*sqrt(5)*(x-2)) + exp(10*sqrt(5)*x))/(1 + exp(20*sqrt(5)))

fig  = figure(figsize=(6 ,4))
ax   = fig.add_subplot(111)

ax.plot(t, uf1, 'r',   lw=2.0, label=r"$u$")
ax.plot(t, uf2, 'k',   lw=2.0, label=r"$\hat{u},\ \tilde{u}$")
ax.plot(x, ue,  'k--', lw=2.0, label=r"$u_{\mathrm{a}}$")

ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$u$')
ax.grid()
leg = ax.legend(loc='upper right')
leg.get_frame().set_alpha(0.0)
#ax.set_xlim([0,0.3])
tight_layout()
show()

---

The actual mistake is in the space on which BC lives:

leftBC  = DirichletBC(M, uD, left)

There must be M, not Q. Workaround to FFC #69 is unimportant here.

---

Oh my, of course. My bad.

---

Function solve should not tolerate this. I've filed an issue at https://bitbucket.org/fenics-project/dolfin/issues/717 and started fixing it.