poisson.py
df534370ac8fdevelopment

poisson.py
View Options

			"""
			FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.
			Test problem is chosen to give an exact solution at all nodes of the mesh.
			-Laplace(u) = f in the unit square
			u = u_D on the boundary
			u_D = 1 + x^2 + 2y^2
			f = -6
			"""

			from __future__ import print_function
			from fenics import *
			import matplotlib.pyplot as plt

			# Create mesh and define function space
			mesh = UnitSquareMesh(8, 8)
			V = FunctionSpace(mesh, 'P', 1)

			# Define boundary condition
			u_D = Expression('1 + x[0]x[0] + 2x[1]*x[1]', degree=2)

			def boundary(x, on_boundary):
			return on_boundary

			bc = DirichletBC(V, u_D, boundary)

			# Define variational problem
			u = TrialFunction(V)
			v = TestFunction(V)
			f = Constant(-6.0)
			a = dot(grad(u), grad(v))*dx
			L = fvdx

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

			# Plot solution and mesh
			plot(u)
			plot(mesh)

			# Save solution to file in VTK format
			vtkfile = File('poisson/solution.pvd')
			vtkfile << u

			# 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
			plt.show()