Hello,
I am trying to make a simple x - y graph to represent heat vs time with my cylindrical heat equation. I am relatively new to fenics/python and would please like some help.
Thanks,
My code and Results:
from future import print_function
from dolfin import *
from mshr import *
from fenics import *
import time
#Making a cylindrical geometry
cylinder = Cylinder(Point(0, 0, -1), Point(0, 0, 1), 0.5, 0.5)
geometry = cylinder
Making Mesh (40 corresponds to the mesh density)
mesh = generate_mesh(geometry, 40)
Save the mesh
File(“cylinder.pvd”) << mesh
V = FunctionSpace(mesh, ‘P’, 1)
T = 2.0 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
Define boundary condition
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, Constant(0), boundary)
Define initial value
u_0 = Expression(‘exp(-apow(x[0], 2) - apow(x[1], 2))’,
degree=2, a=5)
u_n = interpolate(u_0, V)
Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
F = uvdx + dt*dot(grad(u), grad(v))dx - (u_n + dtf)vdx
a, L = lhs(F), rhs(F)
Create VTK file for saving solution
vtkfile = File(‘heat_gaussian/solution.pvd’)
Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(a == L, u, bc)
# Save to file and plot solution
vtkfile << (u, t)
plot(u)
# Update previous solution
u_n.assign(u)
import matplotlib.pyplot as plt
plt.show()