Fluidity PyOP2 identity test case

For reference only: this is a PyOP2 identity test case for running from Fluidity without the convenience layer flop.py.

"""
This is a demo of the use of ffc to generate kernels. It solves the identity
equation on a quadrilateral domain. It requires the fluidity-pyop2 branch of
ffc, which can be obtained with:

bzr branch lp:~grm08/ffc/fluidity-pyop2

This may also depend on development trunk versions of other FEniCS programs.
"""

from pyop2 import op2
from pyop2.ffc_interface import compile_form
from ufl import *

op2.init(backend='sequential')

# Set up finite element identity problem

E = FiniteElement("Lagrange", "triangle", 1)

v = TestFunction(E)
u = TrialFunction(E)
f = Coefficient(E)

a = inner(grad(v),grad(u))*dx + v*u*dx
L = v*f*dx

# Generate code for mass and rhs assembly.

mass_code = compile_form(a, "mass")
rhs_code  = compile_form(L, "rhs")

mass = op2.Kernel(mass_code, "mass_cell_integral_0_0")
rhs  = op2.Kernel(rhs_code,  "rhs_cell_integral_0_0" )

# Set up simulation data structures

t = state.scalar_fields["Tracer"]
c = state.vector_fields["Coordinate"]
valuetype = numpy.float64

#from fluidity_tools import shell
#shell()

nodes = op2.Set(c.node_count, "nodes")
elements = op2.Set(c.element_count, "elements")

elem_node = op2.Map(elements, nodes, 3, c.mesh.ndglno - 1, "elem_node")

sparsity = op2.Sparsity(elem_node, elem_node, 1, "sparsity")
mat = op2.Mat(sparsity, 1, valuetype, "mat")

coords = op2.Dat(nodes, 2, c.val, valuetype, "coords")

x_vals = numpy.zeros(t.node_count)
f = op2.Dat(nodes, 1, t.val, valuetype, "f")
b = op2.Dat(nodes, 1, numpy.zeros(t.node_count), valuetype, "b")
x = op2.Dat(nodes, 1, x_vals, valuetype, "x")

# Assemble and solve

op2.par_loop(mass, elements(3,3),
             mat((elem_node(op2.i(0)), elem_node(op2.i(1))), op2.INC),
             coords(elem_node, op2.READ))

op2.par_loop(rhs, elements,
                     b(elem_node, op2.INC),
                     coords(elem_node, op2.READ),
                     f(elem_node, op2.READ))

op2.solve(mat, b, x)

# Print solution

print "Expected solution: %s" % t.val
print "Computed solution: %s" % x_vals

t.val = x_vals
state.scalar_fields["Tracer"] = t