""" This demo implements a Poisson equations solverbased on the demo "dolfin/demo/pde/poisson/python/demo.py"in Dolfin using Epetra matrices, the AztecOO CG solver and ML AMG preconditioner """__author__ = "Kent-Andre Mardal (kent-and@simula.no)"__date__ = "2008-04-24"__copyright__ = "Copyright (C) 2008 Kent-Andre Mardal"from dolfin import *# Create mesh and finite elementmesh = UnitSquare(20,20)element = FiniteElement("Lagrange", "triangle", 1)# Source termclass Source(Function):    def __init__(self, element, mesh):        Function.__init__(self, element, mesh)    def eval(self, values, x):        dx = x[0] - 0.5        dy = x[1] - 0.5        values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)# Neumann boundary conditionclass Flux(Function):    def __init__(self, element, mesh):        Function.__init__(self, element, mesh)    def eval(self, values, x):        if x[0] > DOLFIN_EPS:            values[0] = 25.0*sin(5.0*DOLFIN_PI*x[1])        else:            values[0] = 0.0# Sub domain for Dirichlet boundary conditionclass DirichletBoundary(SubDomain):    def inside(self, x, on_boundary):        return bool(on_boundary and x[0] < DOLFIN_EPS)# Define variational problemv = TestFunction(element)u = TrialFunction(element)f = Source(element, mesh)g = Flux(element, mesh)a = dot(grad(v), grad(u))*dxL = v*f*dx + v*g*ds# Create backendbackend = EpetraFactory.instance()# Assemble matricesA = assemble(a, mesh, backend=backend)b = assemble(L, mesh, backend=backend) # Define boundary conditionu0 = Function(mesh, 0.0)boundary = DirichletBoundary()bc = DirichletBC(u0, mesh, boundary)bc.apply(A, b, a)# Create solution vector (also used as start vector) x = b.copy()x.zero()#solve(A,x,b, gmres, amg)solve(A,x,b, gmres, jacobi)# plot the solutionU = Function(element, mesh, x)plot(U)interactive()# Save solution to filefile = File("poisson.pvd")file << U