# Examples

Basic program
Basic unstructured mesh generation
Simple nonlinear problem (Newton's Method)
Realistic geometry and unstructured mesh
Simple domain, changing BC's

### Example 1: basic program

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

### Example 2: basic unstructured mesh generation

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

### Example 3: simple nonlinear problem

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 &")

### Example 4: realistic geometry and unstructured mesh

    #!/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()

### Example 5: simple domain, changing BC's

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