How to plot solutions in one variable

Question

Dear all,
I am searching but can't seem to find much information about simulations in one dimension only.
I am trying to make a simple heat transfer model in one dimension.
Please find below the code:

from __future__ import print_function
from fenics import *
import numpy as np
from IPython.display import HTML
import matplotlib.pyplot as plt

T = 10.0 # final timenum_steps = 10 # number of time steps
dt = T / num_steps # time step size`

nx = 20mesh = UnitIntervalMesh(nx)
V = FunctionSpace(mesh, 'P', 1)`

# Define boundary condition
u_D = Expression('t*0.5 + 20.0', degree=2, t=0)
#u_D = Expression('50', degree=1, t=0)

tol = 1e-14
def boundary(x, on_boundary):
return on_boundary and near(x[0], 0.0, tol)

bc = DirichletBC(V, u_D, boundary)

# Define initial value
ic = Expression('20', degree=1)
u_n = interpolate(ic, V)
# u_n = project(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

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
values = []
for n in range(num_steps):

# Update current time
t += dt
u_D.t = t

# Compute solution
solve(a == L, u, bc)

values.append(u.vector().array())
# Update previous solution
u_n.assign(u)

you may see that I have made only minor adaptions on the heat tutorial program.
I am storing the values in a list for plotting with matplotlib, but I think that they are not ordered.
Please any help will be really important since I am a begginer.
(I also know the documention, and Langtangen's book however my mathematical background is not the best yet, this is why I am looking for simple one dimension examples to start)
All the best, Murilo Moreira.

Answer 1

Hi, here's a made up isolated problem to show you what's going on

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

mesh = UnitIntervalMesh(10)

V = FunctionSpace(mesh, 'CG', 2)
f = interpolate(Expression('x[0]*x[0]', degree=2), V)

values = f.vector().array()

plt.figure()
plt.plot(y, linestyle='none', marker='o')  

# As you noticed the coefficients are not ordered according 
# to x coord of dofs. Get the x axis 'right' (align x an values) by e.g.
x = interpolate(Expression('x[0]', degree=1), V).vector().array()
# or
x0 = V.tabulate_dof_coordinates()
assert np.linalg.norm(x-x0) < DOLFIN_EPS  # They are the same thing

plt.figure()
plt.plot(x, values, linestyle='none', marker='o')

# This is better, but what if you wanted a line plot
plt.figure()
plt.plot(x, values)
# --> You see that the plot isn't quite right; the interior dofs of CG2 element are an 
# issue. So finally compute how to allign x according to coordinates
indices = np.argsort(x)
x = x[indices]
values = values[indices] 
plt.figure()
plt.plot(x, values)

# Of course with P1 elements things are simpler and anything following # --> is not 
# necessary 
plt.show()

Comments

Dear MiroK,
Thank you very much!
With your code I could understand how things work!
I am owing you!
Thank you very much for your kind help!
All the best, Murilo.