This example illustrates the use of ellipt2d to solve the Laplace equation on a simple rectangular domain with Dirichlet boundary conditions set only at the top and bottom edges of the domain. It took about ten minutes to write. This example also cleary shows the code structure for solving all other problems with structured meshes. The user may specify a structured mesh and the reg2tri class provides routines to triangulate it. However, one may also lean on a limited number of built in regular meshes (as this example does).
#!/usr/bin/env python
#
import time, os
from math import sqrt
from ellipt2d import ellipt2d
from DirichletBound import DirichletBound
import vector, reg2tri
class demo_Simple:
"""
This demo shows how simple it can be to solve
uncomplicated problems with ellipt2d.
"""
def __init__(self):
# number of nodes
N = 100
# box size
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0
nx1 = int(sqrt(N))
ny1 = int(sqrt(N))
nx1 = max(3, nx1)
ny1 = max(3, ny1)
# Mesh generation.
grid = reg2tri.rect2cross((xmin, ymin, xmax, ymax), nx1, ny1)
# Create Dirichlet boundary condition object
db = DirichletBound()
# Set the boundary conditions
for i in range(0,nx1):
db[i] = 2.0
for i in range( (ny1-1)*(nx1) , (ny1-1)*(nx1)+ nx1 ):
db[i] = 4.0
# Create ellipt2d object.
# Equation is -div f grad v + g v = s
f='1.'; g='1.'; s='0.'
equ = ellipt2d(grid, f, g, s)
# Assemble the stiffness matrix
[amat, rhs] = equ.stiffnessMat()
# Apply the boundary conditions (modify amat and rhs
equ.dirichletB(db,amat,rhs)
v = amat.CGsolve(rhs,rhs,(1.0e-3)/(float(N)), 2*len(rhs) )
# OpenDX plotting
equ.toDX(v,'ellipt2d.dx')
#os.system("dx -image ellipt2d.net")
# Plotting using builtin Tk methods
from Tkinter import Tk, Frame, Button, Canvas, BOTTOM
import tkplot
root = Tk()
frame = Frame(root)
frame.pack()
WIDTH, HEIGHT = 500, 450
button = Button(frame, text="OK", fg="red", command=frame.quit)
button.pack(side=BOTTOM)
canvas = Canvas(bg="white", width=WIDTH, height=HEIGHT)
canvas.pack()
tkplot.tkplot(canvas, grid, v, 0,0,1, WIDTH, HEIGHT)
tkplot.tkplot(canvas, grid, db, 0,0,1, WIDTH, HEIGHT)
root.mainloop()
if __name__ == "__main__":
a = demo_Simple()
This example illustrates the steps neccessary to generate a quality unstructured 2D mesh. Ellipt2d provides an interface with Triangle (see main page), which may be used in standalone mode or within the ellipt2d environment.
#!/usr/bin/env python
from ellipt2d import ellipt2d
from Ireg2tri import Ireg2tri
import vector
import time, os
class demo_Triangulation:
"""
This demonstrates the use of the trinagulation routines..
"""
def __init__(self):
points = []
seglist = []
holelist = []
regionlist = []
points = [(0.0,0.0),(0.2,0.0),(0.6,0.0),(0.8,0.0),(1.0,0.0),(0.0,1.0)]
regionlist.append( (0.2,0.2) )
seglist = [(0,1),(1,2),(2,3),(3,4),(4,5),(5,0)]
initialAreaBound = 0.001
meshgen = Ireg2tri(initialAreaBound)
meshgen.setUpConvexHull(points,seglist,regionlist,holelist)
self.grid = meshgen.triangulate()
self.grid = meshgen.refine(1)
self.grid.plot()
if __name__ == "__main__":
a = demo_Triangulation()
This example takes advantage of the fact that the f, g and s functions defining the elliptic operators
can either be strings
(e.g. 'x**2+y**2') or lists of node values. The latter is useful when solving
quasi-nonlinear problems where these functions depend on the values of the solution on the nodes as
computed in a previous iteration.
Note also that F in the most general
case is a tensor. You can give ((fxx,fxy),(fyx,fyy)), or (fxx, fyy), or fxx for
F in the ellipt2d constructor.
#!/usr/bin/env python
#
# J. Mollis 8/1/2000
from math import sqrt, exp
from ellipt2d import ellipt2d
from DirichletBound import DirichletBound
import reg2tri, vector
import time, os
"""
This demo illustrates the Newton method for solving
non-linear elliptic PDE's using ellipt2d.
"""
N = 100
x0, y0 = 0.0, 0.0
xmax, ymax = 1.0, 1.0
nx1 = int(sqrt(N))
ny1 = int(sqrt(N))
nx1 = max(3, nx1)
ny1 = max(3, ny1)
grid = reg2tri.rect2cross((x0, y0, xmax, ymax), nx1, ny1)
db = DirichletBound()
for i in range(0,nx1):
db[i] = 0.0
for i in range( (ny1-1)*(nx1) , (ny1-1)*(nx1)+ nx1 ):
db[i] = 0.0
for i in range(0,ny1-1):
db[(i)*nx1] = 0.0
db[(i)*nx1 + nx1-1] = 0.0
#--------------
f = '1' # Always equal to one for the duration.
#Initial state ------
g, s = [], []
for i in range(0,len(grid)):
g.append(-1.0)
s.append(-1.0)
v = vector.zeros(len(grid))
ad = 1.0
while ad > 1.0e-6:
equ = ellipt2d(grid, f,g,s)
[amat, s] = equ.stiffnessMat()
equ.dirichletB(db,amat,s)
vn = amat.CGsolve(v, s,(1.0e-3)/(float(N)))
ad = 0.0
for j in range(0,len(v)):
ad = ad + abs(v[j]-vn[j])
g[j] = -exp(-vn[j])
s[j] = -exp(-vn[j]) - vn[j]*exp(-vn[j])
ad = ad/len(v)
v = vn
#--------------
# plot result
from Tkinter import Tk, Frame, Button, Canvas, BOTTOM
import tkplot
root = Tk()
frame = Frame(root)
frame.pack()
WIDTH, HEIGHT = 500, 450
button = Button(frame, text="OK", fg="red", command=frame.quit)
button.pack(side=BOTTOM)
canvas = Canvas(bg="white", width=WIDTH, height=HEIGHT)
canvas.pack()
tkplot.tkplot(canvas, grid, v, 0,0,1, WIDTH, HEIGHT)
tkplot.tkplot(canvas, grid, db, 0,0,1, WIDTH, HEIGHT)
root.mainloop()
# or using openDX
equ.toDX(v, 'ellipt2d.dx' )
#os.system("dx -image ellipt2d.net &")
#!/usr/bin/env python
#
# A. Pletzer, J. Mollis 8/1/2000
from math import pi
from ellipt2d import ellipt2d
from Ireg2tri import Ireg2tri
from DirichletBound import DirichletBound
import vector
import time, os
try:
from geqdsk import geqdsk
except:
raise "this demo requires geqdsk... (not supplied with ellipt2d)"
import superlu
class demo_RealisticIreg:
"""
Unstructured mesh with realistic irregular geometryRe
"""
def __init__(self):
# 1. Equation definition--------------
k = 2.*pi
self.f_funct_str, self.g_funct_str, self.s_funct_str = '1', '%10.4f' % (-k**2), '0'
# 2. Domain grid/boundary defintions----------------
g = geqdsk("./g101544.00185")
pps = g.getPlasmaGeometry()
vacps = g.getLimiterGeometry()
points = []
seglist = []
holelist = []
regionlist = []
dB = DirichletBound()
for i in range(0,len(pps)-1):
points.append(pps[i])
dB[i] = 1.0
for i in range(0,len(vacps)-2):
points.append(vacps[i])
dB[len(pps)-1 + i] = 0.0
for j in range(0,len(pps)-2):
seglist.append((j,j+1))
seglist.append((0,len(pps)-2))
pos = len(pps)-1
for j in range(pos,len(points)-1):
seglist.append((j,j+1))
seglist.append((pos,len(points)-1))
holelist.append( (0.99,0.0) )
regionlist.append( (0.99, 0.0) )
regionlist.append( (0.478,-1.18) )
## Grid creation/boundary set up------------------
initialAreaBound = 0.01
meshgen = Ireg2tri(initialAreaBound)
print 'points = ', points
print 'seglist = ', seglist
meshgen.setUpConvexHull(points,seglist,regionlist,holelist)
meshgen.setUpDirichlet(dB)
self.grid = meshgen.triangulate()
self.grid = meshgen.refine(1)
# update the boundary conditions
meshgen.updateDirichlet(dB)
dB.plot(self.grid)
# 4. Set up the equation.
equ = ellipt2d(self.grid, self.f_funct_str, self.g_funct_str, self.s_funct_str)
# 5. Assemble stiffness matrix and compute source vector
t = time.time()
[amat, s] = equ.stiffnessMat()
print 'Time to assemble stiffness matrix:',(time.time()-t),'Sec.'
# 6. Apply Boundary conditions.
equ.dirichletB(dB,amat,s)
#equ.neumannB(nB,s)
#equ.robinB(rB,amat,s)
# 7. Solve linear system.
to = time.time()
v = superlu.solve(amat,s)
print 'Time to solve Linear System:',(time.time()-to),'Sec.'
# 8. Write DX and UCD formated files.
equ.toUCD(v, 'ellipt2d.inp' )
equ.toDX(v, 'ellipt2d.dx' )
#9. Display results with OpenDX or AVS/Express.
#os.system("dx -image ellipt2d.net")
#10. Plotting using builtin Tk methods
from Tkinter import Tk, Frame, Button, Canvas, BOTTOM
import tkplot
root = Tk()
frame = Frame(root)
frame.pack()
WIDTH, HEIGHT = 500, 450
button = Button(frame, text="OK", fg="red", command=frame.quit)
button.pack(side=BOTTOM)
canvas = Canvas(bg="white", width=WIDTH, height=HEIGHT)
canvas.pack()
tkplot.tkplot(canvas, equ.grid, v, 0,0,1, WIDTH, HEIGHT)
root.mainloop()
if __name__ == "__main__":
a = demo_RealisticIreg()
demo_Rbc.py.
#!/usr/bin/env python
from math import sqrt
from ellipt2d import ellipt2d
from DirichletBound import DirichletBound
from NeumannBound import NeumannBound
import reg2tri
import vector
import time, os
class demo_Simple:
"""
This demo sets some interesting boundary conditions. Interior
Dirichlet BCs and and a few Neumann edges.
"""
def __init__(self):
N = 100
f = '1'
g = '1'
s = '0'
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0
nx1 = int(sqrt(N))
ny1 = int(sqrt(N))
nx1 = max(3, nx1)
ny1 = max(3, ny1)
grid = reg2tri.rect2cross((xmin, ymin, xmax, ymax), nx1, ny1)
db = DirichletBound()
for i in range(0,nx1-3):
db[i] = 2.0
db[nx1-3] = 5.0
db[nx1-2] = 6.0
db[nx1-1] = 7.0
db[25] = 8.0 # inner Dirichlet BC's.
db[55] = 12.0
for i in range( (ny1-1)*(nx1) , (ny1-1)*(nx1)+ nx1 ):
db[i] = 4.0
nb = NeumannBound() # Neumann boundary condition set on one segment.
nb[(nx1,(2)*nx1)] = 3.0
nb[(2*nx1,3*nx1)] = 3.0
equ = ellipt2d(grid, f,g,s)
[amat, s] = equ.stiffnessMat()
equ.dirichletB(db,amat,s)
equ.neumannB(nb,s)
v = amat.CGsolve(s,s,(1.0e-3)/(float(N)), 2*len(s) )
# plot result
from Tkinter import Tk, Frame, Button, Canvas, BOTTOM
import tkplot
root = Tk()
frame = Frame(root)
frame.pack()
WIDTH, HEIGHT = 500, 450
button = Button(frame, text="OK", fg="red", command=frame.quit)
button.pack(side=BOTTOM)
canvas = Canvas(bg="white", width=WIDTH, height=HEIGHT)
canvas.pack()
# plot solution Dirichlet boundary conditions
tkplot.tkplot(canvas, grid, v, 0,0,1, WIDTH, HEIGHT)
tkplot.tkplot(canvas, grid, db, 0,0,1, WIDTH, HEIGHT)
tkplot.tkplot(canvas, grid, nb, 0,0,1, WIDTH, HEIGHT)
root.mainloop()
# or using openDX
equ.toDX(v,'ellipt2d.dx')
#os.system("dx -image ellipt2d.net")
if __name__ == "__main__":
a = demo_Simple()