Hi all, I am a new FEniCS user. I am trying to practice the configuration of FEniCS from the simplest one, one of which is a heat equation. I simply copied and pasted the example code, and changed the input mesh. Then, the results were very weird as follows. When I used the crossed mesh, the error has regionally big values whereas the error was small with the left and right mesh. The code is attached here.

“”"

FEniCS tutorial demo program: Heat equation with Dirichlet conditions.

Test problem is chosen to give an exact solution at all nodes of the mesh.

u’= Laplace(u) + f in the unit square

u = u_D on the boundary

u = u_0 at t = 0

u = 1 + x^2 + alpha*y^2 + \beta*t

f = beta - 2 - 2*alpha

“”"

```
from __future__ import print_function
from fenics import *
import numpy as np
T = 2.0 # final time
num_steps = 10 # number of time steps
dt = T / num_steps # time step size
alpha = 3 # parameter alpha
beta = 1.2 # parameter beta
# Create mesh and define function space
nx = ny = 8
#mesh = UnitSquareMesh(nx, ny,'crossed')
#mesh = UnitSquareMesh(nx, ny,'left')
#mesh = UnitSquareMesh(nx, ny,'right')
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
degree=2, alpha=alpha, beta=beta, t=0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)
# Time-stepping
u = Function(V)
t = 0
vtkfile = File('tutorial.pvd')
for n in range(num_steps):
# Update current time
t += dt
u_D.t = t
# Compute solution
solve(a == L, u, bc)
# Plot solution
plot(u)
vtkfile << (u,t)
# Compute error at vertices
u_e = interpolate(u_D, V)
error = np.abs(u_e.vector()[:] - u.vector()[:]).max()
print('t = %.2f: error = %.3g' % (t, error))
# Update previous solution
u_n.assign(u)
# Hold plot
#interactive()
```

A specific description of code is given here,

This is the error plot over the domain with the crossed mesh