Thank you for your time and patience. This one solves my problem, but I find another problem.

When I add

`mesh2.coordinates().sort(axis=0)`

,

my stiffness matrix become singular. I think this function may have a bug.

The whole code is below

```
from fenics import *
import matplotlib.pyplot as plt
import sympy
def AllExpressions(sigma, cell): # the PDE's exact solution and source term
x, y = sympy.symbols("x[0] x[1]")
u = sympy.exp(x*y) * sympy.sin(pi*x) * sympy.sin(pi*y) # exact solution, satisfies boundary condition
u_x = u.diff(x, 1)
u_y = u.diff(y, 1)
u_xx = u_x.diff(x, 1)
u_yy = u_y.diff(y, 1)
f = -(u_xx + u_yy) \
+ (sigma.values()[0]*u_x + sigma.values()[1]*u_y)
u = Expression(sympy.printing.ccode(u), degree=6, cell=cell)
f = Expression(sympy.printing.ccode(f), degree=4, cell=cell)
return u, f
def solve_coarse_problem():
# ------------------------------------------------- coarse mesh problem -------------------------------------------
p, q = TrialFunction(VH), TestFunction(VH)
AH = assemble(inner(grad(p), grad(q)) * dx + inner(sigma, grad(p)) * q * dx)
bH = assemble(f * q * dx)
bcH.apply(AH)
bcH.apply(bH)
solver = KrylovSolver()
solver.set_operator(AH)
uH = Function(VH)
solver.solve(uH.vector(), bH)
plt.figure()
plot(uH, mode="warp", title="solution on coarse space")
def solve_fine_problem():
# ----------------------------------------------- fine mesh problem -----------------------------------------------
u, v = TrialFunction(Vh), TestFunction(Vh)
Ah = assemble(inner(grad(u), grad(v)) * dx + inner(sigma, grad(u)) * v * dx)
bh = assemble(f * v * dx)
bch.apply(Ah)
bch.apply(bh)
solver = KrylovSolver()
solver.set_operator(Ah)
uh = Function(Vh)
solver.solve(uh.vector(), bh)
plt.figure()
plot(uh, mode="warp", title="solution on fine space")
if __name__ == '__main__':
sigma = Constant((8, 8))
cell = "triangle"
u_exact, f = AllExpressions(sigma, cell)
element = FiniteElement("CG", cell, 1)
coarse = UnitSquareMesh(8, 8)
fine = refine(refine(refine(coarse)))
fine.coordinates().sort(axis=0) # comment this line, I can get the right result!!!
VH = FunctionSpace(coarse, element)
bcH = DirichletBC(VH, Constant(0.0), DomainBoundary())
Vh = FunctionSpace(fine, element)
bch = DirichletBC(Vh, Constant(0.0), DomainBoundary())
solve_coarse_problem()
solve_fine_problem()
plt.show()
```