1
# This demo program solves Poisson's equation
3
# - div grad u(x, y) = f(x, y)
5
# on the unit square with source f given by
7
# f(x, y) = 500*exp(-((x-0.5)^2 + (y-0.5)^2)/0.02)
9
# and boundary conditions given by
14
# using a discontinuous Galerkin formulation (interior penalty method).
16
__author__ = "Kristian B. Oelgaard (k.b.oelgaard@tudelft.nl)"
17
__date__ = "2007-10-02 -- 2007-10-02"
18
__copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
19
__license__ = "GNU LGPL Version 2.1"
23
# Create mesh and finite element
24
mesh = UnitSquare(64, 64)
25
element = FiniteElement("Discontinuous Lagrange", "triangle", 1)
28
class Source(Function):
29
def __init__(self, element, mesh):
30
Function.__init__(self, element, mesh)
31
def eval(self, values, x):
34
values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)
36
# Define variational problem
37
# Test and trial functions
38
v = TestFunction(element)
39
u = TrialFunction(element)
40
f = Source(element, mesh)
42
# Normal component, mesh size and right-hand side
43
n = FacetNormal("triangle", mesh)
44
h = AvgMeshSize("triangle", mesh)
51
a = dot(grad(v), grad(u))*dx \
52
- dot(avg(grad(v)), jump(u, n))*dS \
53
- dot(jump(v, n), avg(grad(u)))*dS \
54
+ alpha/h('+')*dot(jump(v, n), jump(u, n))*dS \
55
- dot(grad(v), mult(u, n))*ds \
56
- dot(mult(v, n), grad(u))*ds \
62
# Solve PDE and plot solution
63
pde = LinearPDE(a, L, mesh)
67
# Save solution to file
68
file = File("poisson.pvd")