I am working on a mixed finite element problem for Poisson’s equation using a BDM1 approximation for the flux. I tested my method on a manufactured solution, with piecewise linear flux, and was getting very large errors in the flux. After testing a few things, I am questioning the result returned by errornorm, particularly how it handles interpolating expressions to Vector DG spaces. Here is what I have coded:
Qe = FiniteElement('BDM',triangle,1)
Ue = FiniteElement('CG',triangle,2)
element = MixedElement([Qe, Ue])
V = FunctionSpace(mesh,element)
# ...
# set up problem
# ...
QU = Function(V)
solve(A,QU.vector(),b)
Q,U = QU.split(True)
class qeExpression(UserExpression):
def eval(self, value, x):
if x[0] <= 0.5 and x[1] <=0.5:
value[0] = -10*x[0]
value[1] = 5.2 - 10.4*x[1]
elif x[0] > 0.5 and x[1] <=0.5:
value[0] = 10.*(x[0]-1)
value[1] = 25*(1 - 2*x[1])
elif x[0] <=0.5 and x[1] > 0.5:
value[0] = -x[0]
value[1] = 1 - 2*x[1]
else:
value[0] = x[0] - 1
value[1] = 1.5 - 3*x[1]
def value_shape(self):
return (2,)
qe = qeExpression(degree = 1)
print(errornorm(qe,Q,'Hdiv')) # equals 2.56
pi_qe = interpolate(qe,FunctionSpace(mesh,'BDM',1))
print(errornorm(pi_qe,Q,'Hdiv')) # equals 2.09e-12
diff = pi_qe.vector().vec().array - Q.vector().vec().array))
print(numpy.max(numpy.abs(diff)) # equals 7.10e-13
W = VectorFunctionSpace(mesh,'DG',4)
pi_qe = interpolate(pi_qe,W)
pi_Q = interpolate(Q,W)
e = Function(W)
e.assign(pi_qe)
e.vector().axpy(-1.0,pi_Q.vector())
print(norm(e,norm_type='Hdiv')) # equals 2.09e-12
pi_qe = interpolate(qe,W) # interpolate expression directly
e = Function(W)
e.assign(pi_qe)
e.vector().axpy(-1.0,pi_Q.vector())
print(norm(e,norm_type='Hdiv')) # equals 2.56 again
Since the values of the degrees of freedom are so close (7.10e-13 at worst) between interpolant and the solution on the BDM space, I am much more inclined to believe that errornorm is incorrect and there is something happening when directly interpolating the expression to a VectorDG space.
Is there a better way to define my piecewise defined expression so that this interpolation works as intended?
I have also tried
qe = qeExpression(element=Qe)
instead of defining the degree, but have the same issue.