Python & Finite Elements 
by Alexander Pletzer 

Listing One
(a)
import reg2tri, vector       
from ellipt2d import *

xmin, xmax = 0., 1.4
ymin, ymax = 0., 1.

N = 101
nx1 = int(sqrt(float(N)*(xmax-xmin)/(ymax-ymin)))
ny1 = int(sqrt(float(N)*(ymax-ymin)/(xmax-xmin)))

(b)
grid = reg2tri.rect2criss((xmin, ymin, xmax, ymax), nx1, ny1)


(c)
db = DirichletBound()               
for i in range(nx1):
    db[i]             = 0.0  # south
    db[nx1*(ny1-1)+i] = 0.0   # north

for i in range(ny1):
    y = grid.y(i*nx1)
    db[i*nx1] = sin(pi*(y-ymin)/(ymax-ymin))   # west

(d)
F='1.'
g, s = '0.', '0.'
equ = ellipt2d(grid, F, g, s)

(e)
[amat, b] = equ.stiffnessMat()

(f)
equ.dirichletB(db,amat,b)

(g)
v0 = b
v = amat.CGsolve(v0, b, 1.e-6, 2*len(b))

(h)
equ.toUCD(v, 'ex1.inp')





1
