← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to sandbox/la/neumann/demo.py

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
Import upstream version 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
sandbox/la/neumann/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
    du/dn(x, y) = 0               
 
12
"""
 
13
 
 
14
from dolfin import *
 
15
 
 
16
# Create mesh and finite element
 
17
mesh = UnitSquare(32, 32)
 
18
element = FiniteElement("Lagrange", "triangle", 1)
 
19
 
 
20
# Source term
 
21
class Source(Function):
 
22
    def __init__(self, element, mesh):
 
23
        Function.__init__(self, element, mesh)
 
24
    def eval(self, values, x):
 
25
        dx = x[0] - 0.5
 
26
        dy = x[1] - 0.5
 
27
        values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)
 
28
 
 
29
# Define variational problem
 
30
v = TestFunction(element)
 
31
u = TrialFunction(element)
 
32
f = Source(element, mesh)
 
33
 
 
34
a = dot(grad(v), grad(u))*dx
 
35
m = v*u*dx
 
36
L = v*f*dx 
 
37
 
 
38
x = Vector() 
 
39
b = assemble(L, mesh) 
 
40
A = assemble(a, mesh) 
 
41
M = assemble(m, mesh) 
 
42
 
 
43
z = b.copy() 
 
44
z.assign(1.0) 
 
45
y = b.copy() 
 
46
y.assign(0.0) 
 
47
 
 
48
#A.mult(b,z)
 
49
c  = z.inner(b)/b.size()
 
50
y.assign(c)
 
51
print "type(b) ", type(b)
 
52
print "type(c) ", type(c)
 
53
b -= y
 
54
 
 
55
A.mult(y, z)
 
56
 
 
57
print "zzzzzzzzzzzzzzz"
 
58
print z.norm()
 
59
 
 
60
 
 
61
 
 
62
 
 
63
solve(A, x, b, gmres, amg)
 
64
 
 
65
b.disp()
 
66
x.disp()
 
67
 
 
68
plot(u)
 
69
u  = Function(element, mesh, x)
 
70
 
 
71
 
 
72
# Save solution to file
 
73
file = File("poisson.pvd")
 
74
file << u
 
75
 
 
76
# Hold plot
 
77
plot(u)
 
78
interactive()