I am trying to plot the gradient of a function I solved to observe the discontinuities in the shape functions. I am getting an error when I run the code.

I made the gradient as dTh and I am trying to plot that. My code is as follows:

from fenics import *

import matplotlib.pyplot as plt

import numpy as np

# Generate mesh

domain = IntervalMesh(1000, 0, 1)

SC = FunctionSpace(domain, ‘CG’, 1)

# Define boundary and apply BC

def L_bc(x, on_boundary):

return near(x[0], 0) and on_boundary

def R_bc(x, on_boundary):

return near(x[0], 1) and on_boundary

bc = [DirichletBC(SC, Constant(0.), L_bc), DirichletBC(SC, Constant(0.), R_bc)]

# Define

#f = Expression(‘exp(x[0]*x[0])’, degree=2)

#C = Expression(‘pow(x[0],3) + 10’, degree=3)

f = Expression(‘exp(-x[0])’, degree=2)

#C = Expression(‘2+x[0]*exp(x[0])*(x[0]-1)’, degree=2)

C = Expression(‘2*(6*x[0]*x[0] - 6*x[0] + 1) + x[0]*x[0] exp(-x[0])(x[0]-1)*(x[0]-1)’, degree=2)

# Define test and trial functions

u = TrialFunction(SC)

v = TestFunction(SC)

# Weak form

#a = dot(grad(w), grad(T))*dx + C*w*T*dx

sol = dot(grad(u), grad(v))*dx + v*f*u*dx

#L = f*w*dx

L = v*C*dx

# Compute solution

Th = Function(SC, name=“displacement_linear_1000”)

solve(sol == L, Th, bc)

sigma = grad(u)

x_plt = np.linspace(0, 1, 1000)

Th_plt = np.fromiter(map(Th, x_plt), dtype=np.double)

plt.plot(x_plt, Th_plt)

plt.show()

########################################################

V = VectorFunctionSpace(domain, “DG”, 1)

dTh = project(grad(Th), V)

dTh_plt = np.fromiter(map(dTh, x_plt), dtype=np.double)

plt.plot(x_plt, dTh_plt)

file = File(“displacement_linear_1000.pvd”)

file << dTh