← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/pde/dg/python/demo.py

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/pde/dg/python/demo.py
 
1
# This demo program solves Poisson's equation
 
2
#
 
3
#     - div grad u(x, y) = f(x, y)
 
4
#
 
5
# on the unit square with source f given by
 
6
#
 
7
#     f(x, y) = 500*exp(-((x-0.5)^2 + (y-0.5)^2)/0.02)
 
8
#
 
9
# and boundary conditions given by
 
10
#
 
11
#     u(x, y)     = 0
 
12
#     du/dn(x, y) = 0
 
13
#
 
14
# using a discontinuous Galerkin formulation (interior penalty method).
 
15
 
 
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"
 
20
 
 
21
from dolfin import *
 
22
 
 
23
# Create mesh and finite element
 
24
mesh = UnitSquare(64, 64)
 
25
element = FiniteElement("Discontinuous Lagrange", "triangle", 1)
 
26
 
 
27
# Source term
 
28
class Source(Function):
 
29
    def __init__(self, element, mesh):
 
30
        Function.__init__(self, element, mesh)
 
31
    def eval(self, values, x):
 
32
        dx = x[0] - 0.5
 
33
        dy = x[1] - 0.5
 
34
        values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)
 
35
 
 
36
# Define variational problem
 
37
# Test and trial functions
 
38
v = TestFunction(element)
 
39
u = TrialFunction(element)
 
40
f = Source(element, mesh)
 
41
 
 
42
# Normal component, mesh size and right-hand side
 
43
n = FacetNormal("triangle", mesh)
 
44
h = AvgMeshSize("triangle", mesh)
 
45
 
 
46
# Parameters
 
47
alpha = 4.0
 
48
gamma = 8.0
 
49
 
 
50
# Bilinear form
 
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 \
 
57
   + gamma/h*v*u*ds
 
58
 
 
59
# Linear form
 
60
L = v*f*dx
 
61
 
 
62
# Solve PDE and plot solution
 
63
pde = LinearPDE(a, L, mesh)
 
64
u = pde.solve()
 
65
plot(u)
 
66
 
 
67
# Save solution to file
 
68
file = File("poisson.pvd")
 
69
file << u