← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

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

  • Committer: Niclas Jansson
  • Date: 2010-10-17 10:21:52 UTC
  • Revision ID: njansson@csc.kth.se-20101017102152-e6u3c9uwxa4z19c8
Minor fixes in configure.ac

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()