-
Notifications
You must be signed in to change notification settings - Fork 2
/
poisson-basic.py
62 lines (52 loc) · 1.58 KB
/
poisson-basic.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
from fenics import *
import math
# Use a built-in function to define a
# structured $M\times N$ mesh. It is possible
# import meshes from external mesh generators,
# but that is outside the scope of this
# introduction.
M = 16; N = 32
mesh = RectangleMesh(Point(1,0),
Point(2,math.pi/2),
M,N)
# Introduce a symbolic variable corresponding
# to coordinates in the mesh.
x = SpatialCoordinate(mesh)
# Coefficient and source function for the
# Poisson equation.
kappa = Constant(1)
f = sin(2*x[1])
# Define a function space of continuous Galerkin
# (CG) functions that are polynomials of degree
# one on each element of the mesh.
V = FunctionSpace(mesh,"CG",1)
# `TrialFunction` and `TestFunction` objects
# are symbolic placeholders used as form
# arguments in the definition of linear and
# bilinear forms.
u = TrialFunction(V)
v = TestFunction(V)
# Define $a$ and $L$ from the weak form of the
# Poisson equation as the UFL objects
# `a` and `L`, using math-like syntax.
a = kappa*dot(grad(u),grad(v))*dx
L = f*v*dx
# Define a boundary condition on the space `V`,
# enforcing that functions are zero on the
# entire boundary.
bc = DirichletBC(V,Constant(0),"on_boundary")
# Solve the variational problem, storing the
# solution in a `Function` `u`, which has an
# associated vector of basis function
# coefficients.
u = Function(V)
solve(a==L,u,bc)
# Interactive visualization with Matplotlib:
from matplotlib import pyplot
plot(u)
pyplot.show()
# Write the solution to a Paraview file,
# with a fixed name to re-use saved
# states.
u.rename("u","u")
File("u.pvd") << u