%matplotlib notebook
from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(4, 4)
mesh
V = FunctionSpace(mesh, 'P', 1)
We define the test functions.
u = TrialFunction(V)
v = TestFunction(V)
Here, the boundary conditions are essential: they are not implicit via the variational formulation.
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
This function gets called with a boolean that says if it's on boundary. So here we return the boolean itself.
def boundary(x, on_boundary):
return on_boundary
Now the boundary conditions.
bc = DirichletBC(V, u_D, boundary)
Now we can define the source term:
# Define variational problem
f = Constant(-6.0)
The variational problem:
a = dot(grad(u), grad(v))*dx
L = f*v*dx
Let's define the unknown.
# Compute solution
u = Function(V)
assemble_system(a, L, bc)
(<dolfin.cpp.la.Matrix; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Matrix > *' at 0x11e1f0a50> >, <dolfin.cpp.la.Vector; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Vector > *' at 0x11e1f0f60> >)
And now, let's solve.
solve(a == L, u, bc)
import matplotlib.pyplot as plt
plt.figure()
plot(u, title='Plot')
plot(mesh)
[<matplotlib.lines.Line2D at 0x1188cdb70>, <matplotlib.lines.Line2D at 0x11e2406a0>]
plt.figure()
plot(mesh)
[<matplotlib.lines.Line2D at 0x11c1cb940>, <matplotlib.lines.Line2D at 0x11c1cbba8>]
plt.figure()
cntour = plot(u, scalarbar=True)
plt.colorbar(cntour)
<matplotlib.colorbar.Colorbar at 0x11cd5e240>
# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')
# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))
# Print errors
print('error_L2 =', error_L2)
print('error_max =', error_max)
# Hold plot
interactive()
error_L2 = 0.032940392293421314 error_max = 1.7763568394e-15
Examining nodal values.
nodal_values_u = u.vector()
nodal_values_u
<dolfin.cpp.la.GenericVector; proxy of <Swig Object of type 'std::shared_ptr< dolfin::GenericVector > *' at 0x11cce94e0> >
array_u = nodal_values_u.array()
array_u
array([ 3. , 2.125 , 3.0625, 1.5 , 2.1875, 3.25 , 1.125 , 1.5625, 2.375 , 3.5625, 1. , 1.1875, 1.75 , 2.6875, 4. , 1.0625, 1.375 , 2.0625, 3.125 , 1.25 , 1.6875, 2.5 , 1.5625, 2.125 , 2. ])
vertex_values_u = u.compute_vertex_values()
vertex_values_u
array([ 1. , 1.0625, 1.25 , 1.5625, 2. , 1.125 , 1.1875, 1.375 , 1.6875, 2.125 , 1.5 , 1.5625, 1.75 , 2.0625, 2.5 , 2.125 , 2.1875, 2.375 , 2.6875, 3.125 , 3. , 3.0625, 3.25 , 3.5625, 4. ])
I can plot the solution over a line by evaluating it on a couple of points.
x = np.linspace(0, 1, num=20)
points = [(_, _) for _ in x]
u_line = np.array([u(point) for point in points])
u_exact = np.array([u_D(point) for point in points])
plt.figure()
plt.plot(x, u_line, '-o')
plt.plot(x, u_exact, '-o')
[<matplotlib.lines.Line2D at 0x11ce02438>]
This is pretty cool. I can just evaluate the solution at the points of interest!