Merge lp:~mapdes/fluidity/firedrake-extrusion into lp:fluidity/firedrake

[Doru Bercea]

Adding etruded mesh handling in Firedrake. Changing the global numbering to handle columns of elements and introducing a new Mesh subclass, ExtrudedMesh, to handle the creation of an extruded mesh as a Lagrange 1 function space.
Tests have been introduced to cover: matrix and RHS assembly, unit cube integration with various function space types and extruded helmoltz.
Functionality to output extruded meshes that can be visualized in Paraview has also been introduced.

[David Ham (david-ham)]

Is there a reason why the new tests use gmsh files for unit squares rather than the intrinsic function?

[Lawrence Mitchell (wence)]

On the Fortran side:

Why do you add a check for lshape%type /= -666 in make_mesh?

In make_global_numbering_new:

extrusion is a logical and should be declared as such:

logical :: extrusion

extrusion = .true. etc...

Why introduce dofs_per_stack? What was wrong with dofs_per. They're the same shape are they not? If not, I don't see how things work.

In topology_dofs: Remove the debugging (commented out) dof_old and halo_dof_old stuff, unless it's supposed to actually do something.

Is it clear why we don't use mesh%nodes = total_dofs in both extruded and unextruded cases? It seems to me like they should be the same.

In the python code:

The magic here is horribly underdocumented. A description of the shape of the key datastructures in comments will be very useful if anyone else is to ever have a chance of fixing any bugs.

is it clear why you insist on building fixed-length lists which you then immediately convert into numpy arrays? numpy.zeros(size, dtype=foo) is much better than numpy.array([0]*size, dtype=foo).

Since the ExtrudedMesh class inherits from Mesh it shouldn't have to reimplement all the superclass methods except in the cases where they are different. Why does it? It's imperative to document what the extrusion kernel does, and so on.

Why does fiat_element_from_ufl_element no longer take a ufl_cell as an argument?

Just because FunctionSpace was underdocumented before, adding new functionality is a perfect time to add a docstring. Also, I think it should be an error to pass in vfamily as well as family if the mesh is not an extruded mesh, just for error checking.

Why is the degree now optional? I thought that family was just a string so it can't contain any extra information. But it seems like the extruded case the family is not a string. This /definitely/ needs documenting. I guess the family is now a cross of function spaces? In any case, there's now no error checking if the degree is not passed in an you've asked for a normal function space. So you need to fix that.

op2.Set instantiation does the right thing if layers is passed as None (it sets layers to 1 internally). So there's no need to do that again before Set instantiation.

Python has a != operator, so can you write foo != bar rather than not (foo == bar)?

Function output in the extruded case looks like it's going to be hideously bad, I note that you comment this. Maybe splat a big warning to the screen for now too?

Test directories should be named firedrake_foo with the test python file named test_firedrake_foo.py (so that we can import more than one test in the same python interpreter). As dham says, don't put mesh files in, use the builtin mesh generation routines.

In the tests, why do you return values as strings? Is it for floating point comparison purposes? If so, don't do that, use numpy.allclose. If not, why?

Why are the tolerances on some existing tests altered?

The helmholtz test comments don't match the code. Also, it doesn't actually /test/ anything. Is it a demo?

review: Needs Fixing

[lp:~mapdes/fluidity/firedrake-extrusion]

4281. By Doru Bercea on 2013-09-30

Fixes to Fortran and Python. New LSW test.

4282. By Doru Bercea on 2013-09-30

Name fix.

4283. By Doru Bercea on 2013-09-30

Changes made to fortran and python. Tests renamed and GMSH removed.

4284. By Doru Bercea on 2013-09-30

Two step test added for extrusion.

4285. By Doru Bercea on 2013-09-30

Tests changed and renamed properly.

4286. By Lawrence Mitchell on 2013-10-01

Merge from firedrake trunk

4287. By Lawrence Mitchell on 2013-10-01

Describe why -666 is a magic number when checking mesh continuity

4288. By Lawrence Mitchell on 2013-10-01

Remove dead commented code

4289. By Lawrence Mitchell on 2013-10-01

Add docstrings for extrusion functions

4290. By Lawrence Mitchell on 2013-10-01

Document make_extruded_coords better

4291. By Lawrence Mitchell on 2013-10-01

Clean up make_extruded_coords a little

Don't build lists only to turn them into numpy arrays.

4292. By Lawrence Mitchell on 2013-10-01

Document ExtrudedMesh better

4293. By Lawrence Mitchell on 2013-10-01

Use public API to get a numpy array data pointer

4294. By Lawrence Mitchell on 2013-10-01

Remove duplicate layers property from ExtrudedMesh

4295. By Lawrence Mitchell on 2013-10-01

Document FunctionSpace better

Describe the different methods of getting an extruded function space
as well as a normal function space, and raise errors if the arguments
mismatch in some way.

4296. By Lawrence Mitchell on 2013-10-01

Put back insertion of tensor into mat_cache

4297. By Lawrence Mitchell on 2013-10-01

Whitespace fixup

4298. By Lawrence Mitchell on 2013-10-01

Extruded LSW test doesn't test anything so mark it as special

4299. By Lawrence Mitchell on 2013-10-01

Rename more tests appropriately

4300. By Lawrence Mitchell on 2013-10-01

And final test rename fix

[Lawrence Mitchell (wence)]

I'm now happy with this, but it would be nice to have another pair of eyes.

However, it can't land until all the other bits are in place, which will necessitate kicking Andrew into shape. In particular, all the extrusion tests fail if you have the wrong ffc/ufl/fiat setup, and with the right setup for that, the facet integrals test fails (because of a change Andrew has made in ffc).

review: Approve

[Florian Rathgeber (florian-rathgeber)]

I think this is looking good now.

review: Approve

[lp:~mapdes/fluidity/firedrake-extrusion]

4301. By Florian Rathgeber on 2013-10-03

Merge Firedrake trunk

4302. By Florian Rathgeber on 2013-10-03

Fix typo in ufl.OuterProductElement

4303. By Florian Rathgeber on 2013-10-03

Mark extrusion firedrake tests as such

=== modified file 'femtools/Fields_Allocates.F90'
--- femtools/Fields_Allocates.F90	2013-09-18 14:33:36 +0000
+++ femtools/Fields_Allocates.F90	2013-10-03 16:15:40 +0000
@@ -989,7 +989,7 @@
 
   end function wrap_tensor_field
 
-  function make_mesh (model, shape, continuity, name, with_faces) &
+  function make_mesh (model, shape, continuity, name, with_faces, dofs_per_column) &
        result (mesh)
     !!< Produce a mesh based on an old mesh but with a different shape and/or continuity.
     type(mesh_type) :: mesh
@@ -999,6 +999,7 @@
     integer, intent(in), optional :: continuity
     character(len=*), intent(in), optional :: name
     logical, intent(in), optional :: with_faces
+    integer, dimension(0:), intent(in), optional :: dofs_per_column
 
     integer, dimension(:), allocatable :: ndglno
     real, dimension(:), pointer :: val
@@ -1043,8 +1044,10 @@
     mesh%periodic=model%periodic
 
     ! You can't have a CG degree 0 mesh!
+    ! Fluidity doesn't know about outer product elements (lshape%type
+    ! == -666) so complain.
     if(mesh%shape%degree==0.and.mesh%continuity>=0.and.mesh%shape&
-         &%type/=ELEMENT_TRACE) then
+         &%type/=ELEMENT_TRACE.and.lshape%type/=-666) then
        FLExit("For a P0 mesh, the 'mesh_continuity' must be Discontinuous.")
     end if
 
@@ -1054,7 +1057,8 @@
     end if
 
     mesh%element_classes = model%element_classes
-    call make_global_numbering_new(mesh)
+
+    call make_global_numbering_new(mesh, dofs_per_column=dofs_per_column)
 
     if (associated(model%element_halos)) then
        assert(element_halo_count(model) > 0)

=== modified file 'femtools/Global_Numbering.F90'
--- femtools/Global_Numbering.F90	2013-08-06 12:52:44 +0000
+++ femtools/Global_Numbering.F90	2013-10-03 16:15:40 +0000
@@ -56,10 +56,11 @@
 
 contains
 
-  subroutine make_global_numbering_new(mesh)
+  subroutine make_global_numbering_new(mesh, dofs_per_column)
     ! Construct a new global numbering for mesh using the ordering from
     !  topology.
     type(mesh_type), intent(inout), target :: mesh
+    integer, intent(in), dimension(0:), optional :: dofs_per_column
 
     type(csr_sparsity) :: facet_list, edge_list
     type(csr_matrix) :: facet_numbers, edge_numbers
@@ -70,7 +71,8 @@
     ! The universal identifier. Required for parallel.
     type(scalar_field), pointer :: uid
     integer :: max_vertices, total_dofs, start, finish, d
-    logical :: have_facets, have_halos
+    logical :: have_facets, have_halos, extrusion
+    integer, dimension(0:mesh_dim(mesh)) :: dofs_per_stack
 
     type(integer_set), dimension(:,:), allocatable :: entity_send_targets
     integer, dimension(:), allocatable :: entity_owner
@@ -127,6 +129,15 @@
     ! Extract halo information for each topological entity.
     call create_topology_halos
 
+    ! In the case of extruded meshes change the numbering
+    if (present(dofs_per_column)) then
+        extrusion=.True.
+        dofs_per_stack=dofs_per_column
+    else
+        extrusion=.False.
+        dofs_per_stack=dofs_per
+    end if
+
     ! Save core/non-core/l1/l2 distinction in mesh node_classes
     if (have_halos) then
        
@@ -137,18 +148,18 @@
           finish = finish + entity_counts(d)
           ! Count core nodes (primary sort key -2)
           mesh%node_classes(CORE_ENTITY) = mesh%node_classes(CORE_ENTITY) + &
-               & count(entity_sort_list(start:finish, 0) == -2)*dofs_per(d)
+               & count(entity_sort_list(start:finish, 0) == -2)*dofs_per_stack(d)
           ! Count non-core nodes (primary sort key -1)
           mesh%node_classes(NON_CORE_ENTITY) = &
                & mesh%node_classes(NON_CORE_ENTITY) + &
-               & count(entity_sort_list(start:finish, 0) == -1)*dofs_per(d)
+               & count(entity_sort_list(start:finish, 0) == -1)*dofs_per_stack(d)
           ! L1 halo (primary sort key \in [0, COMM_SIZE]) 
           mesh%node_classes(EXEC_HALO_ENTITY) = mesh%node_classes(EXEC_HALO_ENTITY) + &
                & count(entity_sort_list(start:finish, 0) >= 0 .and. &
-               & entity_sort_list(start:finish, 0) <= getnprocs())*dofs_per(d)
+               & entity_sort_list(start:finish, 0) <= getnprocs())*dofs_per_stack(d)
           ! L2 halo (primary sort key \in [COMM_SIZE + 1, \infty))
           mesh%node_classes(NON_EXEC_HALO_ENTITY) = mesh%node_classes(NON_EXEC_HALO_ENTITY) + &
-               & count(entity_sort_list(start:finish, 0) > getnprocs())*dofs_per(d)
+               & count(entity_sort_list(start:finish, 0) > getnprocs())*dofs_per_stack(d)
        end do
     end if
 
@@ -514,10 +525,10 @@
       halo_dof=0
       do i=1,size(visit_order)
          entity_dof_starts(visit_order(i))=dof
-         dof=dof+dofs_per(entity_dim(visit_order(i)))
+         dof=dof+dofs_per_stack(entity_dim(visit_order(i)))
 
          halo_entity_dof_starts(halo_visit_order(i))=halo_dof
-         halo_dof=halo_dof+dofs_per(entity_dim(halo_visit_order(i)))         
+         halo_dof=halo_dof+dofs_per_stack(entity_dim(halo_visit_order(i)))
       end do
       
       assert(halo_dof == dof)
@@ -582,8 +593,8 @@
                ele_dofs(element%entity2dofs(d,e)%dofs)=&
                     entity_dofs(d, entity, dofs_per)
                
-               dof_to_halo_dof(entity_dofs(d, entity, dofs_per)) &
-                    = halo_entity_dofs(d, entity, dofs_per)
+               dof_to_halo_dof(entity_dofs(d, entity, dofs_per_stack)) &
+                    = halo_entity_dofs(d, entity, dofs_per_stack)
 
             end do
          end do
@@ -594,12 +605,15 @@
       
       assert(all(dof_to_halo_dof>0))
 
-      if (size(mesh%ndglno)>0) then
-         mesh%nodes=maxval(mesh%ndglno)
+      if (extrusion) then
+        mesh%nodes=total_dofs
       else
-         mesh%nodes=0
+        if (size(mesh%ndglno)>0) then
+          mesh%nodes=maxval(mesh%ndglno)
+        else
+          mesh%nodes=0
+        end if
       end if
-         
     end subroutine populate_ndglno
 
     subroutine create_halos(mesh, thalos)
@@ -630,7 +644,7 @@
          
          if (entity_owner(entity)==rank) then
             ! One of ours
-            local_dofs=local_dofs+dofs_per(entity_dim(entity))
+            local_dofs=local_dofs+dofs_per_stack(entity_dim(entity))
             
             do h=1,size(mesh%halos)
                do k=1,key_count(entity_send_targets(entity,h))
@@ -638,7 +652,7 @@
                   if (fetch(entity_send_targets(entity,h),k)==rank) cycle
 
                   call insert(sends(h,fetch(entity_send_targets(entity,h),k)),&
-                       & entity_dofs(entity_dim(entity),entity, dofs_per))
+                       & entity_dofs(entity_dim(entity),entity, dofs_per_stack))
                   
                end do
             end do
@@ -651,7 +665,7 @@
             do h=size(mesh%halos),entity_receive_level(entity),-1
                
                call insert(receives(h,entity_owner(entity)),&
-                    entity_dofs(entity_dim(entity),entity, dofs_per))
+                    entity_dofs(entity_dim(entity),entity, dofs_per_stack))
 
             end do
          end if

=== modified file 'femtools/Python_Interface_f.F90'
--- femtools/Python_Interface_f.F90	2013-09-19 12:59:14 +0000
+++ femtools/Python_Interface_f.F90	2013-10-03 16:15:40 +0000
@@ -197,17 +197,52 @@
     end subroutine populate_facet_maps
 
   end function function_space_f
-
+  
+  function extruded_mesh_f(mesh_ptr, fiat_element, dofs_per_column_ptr) bind (c)
+    type(c_ptr), value, intent(in) :: mesh_ptr
+    type(element_t) :: fiat_element
+    type(function_space_t) :: extruded_mesh_f
+
+    type(mesh_type), pointer :: mesh, fluidity_function_space
+    type(element_type) :: element
+    type(c_ptr), intent(in), value :: dofs_per_column_ptr
+    integer, dimension(:), pointer:: dofs_per_column
+    
+    call c_f_pointer(mesh_ptr, mesh)
+    
+    allocate(fluidity_function_space)
+    
+    element = make_element_shape_from_fiat(fiat_element)
+
+    call c_f_pointer(dofs_per_column_ptr, dofs_per_column, [element%dim+1])
+
+    fluidity_function_space = make_mesh(mesh, element, dofs_per_column=dofs_per_column, with_faces=.false.)
+    
+    extruded_mesh_f%fluidity_mesh = c_loc(fluidity_function_space)
+
+    extruded_mesh_f%element_dof_list = fluidity_function_space%ndglno_c
+
+    extruded_mesh_f%element_count = fluidity_function_space%elements
+    
+    extruded_mesh_f%dof_count = fluidity_function_space%nodes
+
+    extruded_mesh_f%dof_classes = c_loc(fluidity_function_space%node_classes)    
+
+  end function extruded_mesh_f
+  
+  
   function halo_f(mesh_ptr, halo_entity_type) bind(c)
     type(c_ptr), value, intent(in) :: mesh_ptr
+    
     integer(c_int), value, intent(in) :: halo_entity_type
     type(halo_t) :: halo_f
     type(mesh_type), pointer :: mesh
     type(halo_type), pointer :: halo => null()
     integer :: i
     integer, dimension(:), pointer :: tmp
+    
     call c_f_pointer(mesh_ptr, mesh)
-
+    
     select case (halo_entity_type)
     case (VERTEX)
        if (associated(mesh%halos)) then
@@ -412,7 +447,7 @@
     integer :: d, e, f, dof
     integer, allocatable :: count(:,:)
     integer, pointer :: dofs_per(:), entity_dofs(:,:)
-    
+
     element%cell => find_cell(fiat_element%dimension, fiat_element%vertices)
 
     element%dim = fiat_element%dimension

=== modified file 'python/firedrake/core_types.pyx'
--- python/firedrake/core_types.pyx	2013-09-19 14:10:55 +0000
+++ python/firedrake/core_types.pyx	2013-10-03 16:15:40 +0000
@@ -67,10 +67,149 @@
     if not initialised:
         op2.init()
         initialised = True
-
-def fiat_from_ufl_element(ufl_element, ufl_cell):
-    return FIAT.supported_elements[ufl_element.family()]\
-            (_FIAT_cells[ufl_cell.cellname()](), ufl_element.degree())
+            
+def fiat_from_ufl_element(ufl_element):
+    if isinstance(ufl_element, ufl.OuterProductElement):
+        a = FIAT.supported_elements[ufl_element._A.family()]\
+            (_FIAT_cells[ufl_element._A.cell().cellname()](), ufl_element._A.degree())
+
+        b = FIAT.supported_elements[ufl_element._B.family()]\
+            (_FIAT_cells[ufl_element._B.cell().cellname()](), ufl_element._B.degree())
+
+        return FIAT.TensorFiniteElement(a, b)
+    else:
+        return FIAT.supported_elements[ufl_element.family()]\
+            (_FIAT_cells[ufl_element.cell().cellname()](), ufl_element.degree())
+
+# Functions related to the extruded case
+def compute_extruded_dofs(fiat_element, flat_dofs, layers):
+    """Compute the number of dofs in a column"""
+    size = len(flat_dofs)
+    dofs_per_column = np.zeros(size, np.int32)
+    for i in range(size):
+        for j in range(2): #2 is due to the process of extrusion
+            dofs_per_column[i] += (layers - j) * len(fiat_element.entity_dofs()[(i,j)][0])
+    return dofs_per_column
+    
+def compute_vertical_offsets(ent_dofs, flat_dofs):
+    """Compute the offset between corresponding dofs in layers.
+
+    offsets[i] is the offset from the bottom of the stack to the
+    corresponding dof in the ith layer.
+    """
+    size = len(flat_dofs)
+    offsets_per_vertical = np.zeros(size, np.int32)
+    for i in range(size):
+        if len(flat_dofs[i][0]) > 0:
+            offsets_per_vertical[i] = len(flat_dofs[i][0]) - len(ent_dofs[(i,0)][0])
+    return offsets_per_vertical
+    
+def compute_offset(ent_dofs, flat_dofs, total_dofs):
+    """Compute extruded offsets for flattened element.
+
+    offsets[i] is the number of dofs in the vertical for the ith
+    column of flattened mesh entities."""
+    size = len(flat_dofs)
+    res = np.zeros(total_dofs, np.int32)
+    vert_dofs = compute_vertical_offsets(ent_dofs, flat_dofs)
+    for i in range(size):
+        elems = len(flat_dofs[i])
+        dofs_per_elem = len(flat_dofs[i][0])
+        for j in range(elems):
+            for k in range(dofs_per_elem):
+                res[flat_dofs[i][j][k]] = vert_dofs[i]
+    return res
+    
+def total_num_dofs(flat_dofs):
+    """Compute the total number of degrees of freedom in the extruded mesh"""
+    size = len(flat_dofs)
+    total = 0
+    for i in range(size):
+        total += len(flat_dofs[i]) * len(flat_dofs[i][0])
+    return total
+    
+def make_flat_fiat_element(ufl_cell_element, ufl_cell, flattened_entity_dofs):
+    """Create a modified FIAT-style element.
+    Transform object from 3D-Extruded to 2D-flattened FIAT-style object."""
+    # Create base element
+    base_element = fiat_from_ufl_element(ufl_cell_element)
+
+    # Alter base element
+    base_element.dual.entity_ids = flattened_entity_dofs
+    base_element.poly_set.num_members = total_num_dofs(flattened_entity_dofs)
+
+    return base_element
+    
+def make_extruded_coords(mesh, kernel_str, layers):
+    """Given a kernel, use it to generate the extruded coordinates.
+
+    :arg mesh: the 2d mesh to extrude
+    :arg layers: the number of layers in the extruded mesh
+    :arg kernel_str: a C function which produces the z value
+    of the coordinates.  Its calling signature is:
+
+        void extrusion_kernel(double *extruded_coords[],
+                              double *two_d_coords[],
+                              int *layer_number[])
+
+    So for example to build an evenly-spaced extruded mesh with eleven
+    layers and height 1 in the vertical you would write:
+
+       void extrusion_kernel(double *extruded_coords[],
+                             double *two_d_coords[],
+                             int *layer_number[]) {
+           extruded_coords[0][0] = two_d_coords[0][0]; // X
+           extruded_coords[0][1] = two_d_coords[0][1]; // Y
+           extruded_coords[0][2] = 0.1 * layer_number[0][0]; // Z
+       }
+    """
+    coords_dim = len(mesh._coordinates[0])
+    coords_xtr_dim = 3 #dimension
+    kernel = op2.Kernel(kernel_str, "extrusion_kernel")
+    # BIG TRICK HERE:
+    # We need the +1 in order to include the entire column of vertices.
+    # Extrusion is meant to iterate over the 3D cells which are layer - 1 in number.
+    # The +1 correction helps in the case of iteration over vertices which need
+    # one extra layer.
+    iterset = op2.Set(mesh.num_vertices(), "verts1", layers=(layers+1))
+    vnodes = op2.DataSet(iterset, coords_dim)
+    lnodes = op2.DataSet(iterset, 1)
+    nodes_xtr = op2.Set(mesh.num_vertices()*layers, "verts_xtr")
+    d_nodes_xtr = op2.DataSet(nodes_xtr, coords_xtr_dim)
+    d_lnodes_xtr = op2.DataSet(nodes_xtr, 1)
+
+    # Create an op2.Dat with the base mesh coordinates
+    coords_vec = mesh._coordinates.flatten()
+    coords = op2.Dat(vnodes, coords_vec, np.float64, "dat1")
+
+    # Create an op2.Dat with slots for the extruded coordinates
+    coords_new = np.empty(layers * mesh.num_vertices() * coords_xtr_dim, dtype=np.float64)
+    coords_xtr = op2.Dat(d_nodes_xtr, coords_new, np.float64, "dat_xtr")
+
+    # Creat an op2.Dat to hold the layer number
+    layer_vec = np.tile(np.arange(0, layers), mesh.num_vertices())
+    layer = op2.Dat(d_lnodes_xtr, layer_vec, np.int32, "dat_layer")
+
+    # Map a map for the bottom of the mesh.
+    vertex_to_coords = [ i for i in range(0, mesh.num_vertices()) ]
+    v2coords_offset = np.zeros(1, np.int32)
+    map_2d = op2.Map(iterset, iterset, 1, vertex_to_coords, "v2coords", v2coords_offset)
+
+    # Create Map for extruded vertices
+    vertex_to_xtr_coords = [ layers * i for i in range(0, mesh.num_vertices()) ]
+    v2xtr_coords_offset = np.array([1], np.int32)
+    map_xtr = op2.Map(iterset, nodes_xtr, 1, vertex_to_xtr_coords, "v2xtr_coords", v2xtr_coords_offset)
+
+    # Create Map for layer number
+    v2xtr_layer_offset = np.array([1], np.int32)
+    layer_xtr = op2.Map(iterset, nodes_xtr, 1, vertex_to_xtr_coords, "v2xtr_layer", v2xtr_layer_offset)
+
+    op2.par_loop(kernel, iterset,
+                 coords_xtr(op2.INC, map_xtr),
+                 coords(op2.READ, map_2d),
+                 layer(op2.READ, layer_xtr))
+
+    return coords_xtr.data
 
 # C utility functions, not seen in python module
 cdef void function_space_destructor(object capsule):
@@ -88,22 +227,27 @@
     descriptions.
     """
     cdef ft.element_t element
-    topology = fiat_element.ref_el.topology
-    element.dimension = len(topology) - 1
-    element.vertices = len(fiat_element.ref_el.vertices)
-    element.ndof = fiat_element.space_dimension()
-    element.degree = fiat_element.degree()
+    if isinstance(fiat_element, FIAT.tensor_finite_element.TensorFiniteElement):
+        # Use the reference element of the horizontal element.
+        f_element = fiat_element.flattened_element()
+    else:
+        f_element = fiat_element
+        
+    element.dimension = len(f_element.ref_el.topology) - 1
+    element.vertices = len(f_element.ref_el.vertices)
+    element.ndof = f_element.space_dimension()
+    element.degree = f_element.degree()
 
     # Need to free these when you're done
     element.dofs_per = <int *>malloc((element.dimension + 1) * sizeof(int))
     element.entity_dofs = <int *>malloc(3 * element.ndof * sizeof(int))
 
-    entity_dofs = fiat_element.entity_dofs()
+    entity_dofs = f_element.entity_dofs()
     cdef int d
     cdef int e
     cdef int dof
     cdef int dof_pos = 0
-    for d, dim_entities in entity_dofs.iteritems():
+    for d, dim_entities in entity_dofs.iteritems():     
         element.dofs_per[d] = 0
         for e, this_entity_dofs in dim_entities.iteritems():
             # This overwrites the same data a few times, but all of
@@ -115,6 +259,7 @@
                 element.entity_dofs[dof_pos+1] = e
                 element.entity_dofs[dof_pos+2] = dof
                 dof_pos += 3
+    
     return element
 
 class _Facets(object):
@@ -160,6 +305,8 @@
     def __init__(self, *args):
 
         _init()
+        
+        self._layers = 1
 
         self.cell_halo = None
         self.vertex_halo = None
@@ -257,6 +404,10 @@
         for measure in [ufl.dx, ufl.ds, ufl.dS]:
             measure._domain_data = self._coordinate_field
 
+    @property
+    def layers(self):
+        return self._layers
+    
     def cells(self):
         return self._cells
 
@@ -299,6 +450,89 @@
         '''Currently a no-op for flop.py compatibility.'''
         pass
 
+class ExtrudedMesh(Mesh):
+    """Build an extruded mesh from a 2D input mesh
+
+    :arg mesh: 2D unstructured mesh
+    :arg layers: number of structured layers in the "vertical"
+    direction
+    :arg kernel: C kernel to produce 3D coordinates for the extruded
+    mesh see :func:`make_extruded_coords` for more details."""
+    def __init__(self, mesh, layers, kernel):
+        self._old_mesh = mesh
+        self._layers = layers
+        self._cells = mesh._cells
+        self.cell_halo = mesh.cell_halo
+        self._entities = mesh._entities
+        self.parent = mesh.parent
+        self.uid = mesh.uid
+        self.region_ids = mesh.region_ids
+        self.cell_classes = mesh.cell_classes
+        self._coordinates = mesh._coordinates
+        self.name = mesh.name
+
+        self.ufl_cell_element = ufl.FiniteElement("Lagrange",
+                                               domain = mesh._ufl_cell,
+                                               degree = 1)
+        self.ufl_interval_element = ufl.FiniteElement("Lagrange",
+                                               domain = ufl.Cell("interval",1),
+                                               degree = 1)
+
+        self.fiat_base_element = fiat_from_ufl_element(self.ufl_cell_element)
+        self.fiat_vert_element = fiat_from_ufl_element(self.ufl_interval_element)
+        
+        fiat_element = FIAT.tensor_finite_element.TensorFiniteElement(self.fiat_base_element, self.fiat_vert_element)
+        
+        self._ufl_cell = ufl.OuterProductCell(mesh._ufl_cell, ufl.Cell("interval",1))
+        
+        flat_temp = fiat_element.flattened_element()
+        
+        # Calculated dofs_per_column from flattened_element and layers.
+        # The mirrored elements have to be counted only once.
+        # Then multiply by layers and layers - 1 accordingly.
+        self.dofs_per_column = compute_extruded_dofs(fiat_element, flat_temp.entity_dofs(), layers)
+        
+        #Compute Coordinates of the extruded mesh
+        self._coordinates = make_extruded_coords(mesh, kernel, layers)
+
+        # Now we need to produce the extruded mesh using 
+        # techqniues employed when computing the 
+        # function space.
+        cdef ft.element_t element_f = as_element(fiat_element)
+        cdef void *fluidity_mesh = py.PyCapsule_GetPointer(mesh._fluidity_mesh, _mesh_name)
+        if fluidity_mesh == NULL:
+            raise RuntimeError("Didn't find fluidity mesh pointer in mesh %s" % mesh)
+
+        cdef int _ptr = self.dofs_per_column.ctypes.data
+        cdef int *dofs_per_column = <int *><void *>_ptr
+        
+        extruded_mesh = ft.extruded_mesh_f(fluidity_mesh, &element_f, dofs_per_column)
+        
+        #Assign the newly computed extruded mesh as the fluidity mesh
+        self._fluidity_mesh = py.PyCapsule_New(extruded_mesh.fluidity_mesh,
+                                                         _mesh_name,
+                                                         NULL)
+
+        self._coordinate_fs = VectorFunctionSpace(self, "Lagrange", 1)
+
+        self._coordinate_field = Function(self._coordinate_fs,
+                                          val = self._coordinates)
+
+
+        # Set the domain_data on all the default measures to this coordinate field.
+        for measure in [ufl.dx, ufl.ds, ufl.dS]:
+            measure._domain_data = self._coordinate_field
+
+    @utils.cached_property
+    def cell_set(self):
+        if self.cell_halo:
+            size = self.cell_classes
+            halo = self.cell_halo.op2_halo
+        else:
+            size = self.num_cells()
+            halo = None
+        return self.parent.cell_set if self.parent else \
+            op2.Set(size, "%s_elements" % self.name, halo=halo, layers=self._layers)
 
 class Halo(object):
     """Fluidity Halo type"""
@@ -428,24 +662,116 @@
 
 
 class FunctionSpace(object):
-    def __init__(self, mesh, family, degree, name = None):
-
-        self._ufl_element = ufl.FiniteElement(family,
-                                              domain = mesh._ufl_cell,
-                                              degree = degree)
-
-        self.fiat_element = fiat_from_ufl_element(self._ufl_element,
-                                                  mesh._ufl_cell)
-
+    """Create a function space
+
+    :arg mesh: mesh to build the function space on
+    :arg family: string describing function space family, or a
+    :class:`ufl.OuterProductElement`
+    :arg degree: degree of the function space
+    :arg name: (optional) name of the function space
+    :arg vfamily: family of function space in vertical dimension
+    (:class:`ExtrudedMesh`\es only)
+    :arg vdegree: degree of function space in vertical dimension
+    (:class:`ExtrudedMesh`\es only)
+
+    If the mesh is an :class:`ExtrudedMesh`, and the `family` argument
+    is a :class:`ufl.OuterProductElement`, `degree`, `vfamily` and
+    `vdegree` are ignored, since the `family` provides all necessary
+    information, otherwise a :class:`ufl.OuterProductElement` is built
+    from the (`family`, `degree`) and (`vfamily`, `vdegree`) pair.  If
+    the `vfamily` and `vdegree` are not provided, the vertical element
+    will be the same as the provided (`family`, `degree`) pair.
+
+    If the mesh is not an :class:`ExtrudedMesh`, the `family` must be
+    a string describing the finite element family to use, and the
+    `degree` must be provided, `vfamily` and `vdegree` are ignored in
+    this case.
+    """
+    def __init__(self, mesh, family, degree=None, name=None, vfamily=None, vdegree=None):       
+        if isinstance(mesh, ExtrudedMesh):
+            # The function space must adapt to the extrusion
+            # accordingly.
+            # The bottom layer maps  will come from lement_dof_list
+            # dof_count is the total number of dofs in the extruded mesh
+
+            if isinstance(family, ufl.OuterProductElement):
+
+                a = family._A
+                b = family._B
+
+                la = ufl.FiniteElement(a.family(), 
+                                       domain=mesh._old_mesh._ufl_cell, 
+                                       degree=a.degree())
+
+                lb = ufl.FiniteElement(b.family(), 
+                                       domain=ufl.Cell("interval",1), 
+                                       degree=b.degree())
+                                              
+            else:
+                la = ufl.FiniteElement(family, 
+                                       domain=mesh._old_mesh._ufl_cell, 
+                                       degree=degree)
+                # FIAT version of the extrusion
+                if vfamily is None or vdegree is None:
+                    lb = ufl.FiniteElement(family, 
+                                           domain=ufl.Cell("interval",1), 
+                                           degree=degree)
+                else:
+                    lb = ufl.FiniteElement(vfamily, 
+                                           domain=ufl.Cell("interval",1), 
+                                           degree=vdegree)
+                                              
+            # Create the Function Space element
+            self._ufl_element = ufl.OuterProductElement(la, lb)
+
+            # Compute the FIAT version of the UFL element above
+            self.fiat_element = fiat_from_ufl_element(self._ufl_element)
+
+            # Get the flattened version of the 3D FIAT element
+            flat_temp = self.fiat_element.flattened_element()
+
+            # Compute the dofs per column
+            self.dofs_per_column = compute_extruded_dofs(self.fiat_element,
+                                                         flat_temp.entity_dofs(),
+                                                         mesh._layers)
+
+            # Compute the offset for the extrusion process
+            self.offset = compute_offset(self.fiat_element.entity_dofs(),
+                                         flat_temp.entity_dofs(),
+                                         self.fiat_element.space_dimension())            
+        else:
+            if isinstance(family, ufl.OuterProductElement):
+                raise RuntimeError("You can't build an extruded element on an unextruded Mesh")
+            if degree is None:
+                raise RuntimeError("The function space requires a degree")
+
+            self.offset = None
+            self.dofs_per_column = np.zeros(1, np.int32)
+            self.extruded = False
+            
+            self._ufl_element = ufl.FiniteElement(family,
+                                                  domain=mesh._ufl_cell,
+                                                  degree=degree)
+
+            self.fiat_element = fiat_from_ufl_element(self._ufl_element)
+   
+        # Create the extruded function space
         cdef ft.element_t element_f = as_element(self.fiat_element)
+ 
         cdef void *fluidity_mesh = py.PyCapsule_GetPointer(mesh._fluidity_mesh, _mesh_name)
         if fluidity_mesh == NULL:
             raise RuntimeError("Didn't find fluidity mesh pointer in mesh %s" % mesh)
-        
-        function_space = ft.function_space_f(fluidity_mesh, &element_f)
+       
+        cdef int *dofs_per_column = <int *>np.PyArray_DATA(self.dofs_per_column)
+     
+        if isinstance(mesh, ExtrudedMesh):
+            function_space = ft.extruded_mesh_f(fluidity_mesh, &element_f, dofs_per_column)
+        else:      
+            function_space = ft.function_space_f(fluidity_mesh, &element_f)
 
         free(element_f.dofs_per)
         free(element_f.entity_dofs)
+        
         self._fluidity_function_space = py.PyCapsule_New(function_space.fluidity_mesh,
                                                          _function_space_name,
                                                          &function_space_destructor)
@@ -464,13 +790,14 @@
 
         self._dim = 1
 
-        self.interior_facet_node_list = \
-            np.array(<int[:self._mesh.interior_facets.count,:2*element_f.ndof]>
-                      function_space.interior_facet_node_list)
+        if not isinstance(self._mesh, ExtrudedMesh):
+            self.interior_facet_node_list = \
+                np.array(<int[:self._mesh.interior_facets.count,:2*element_f.ndof]>
+                         function_space.interior_facet_node_list)
     
-        self.exterior_facet_node_list = \
-            np.array(<int[:self._mesh.exterior_facets.count,:element_f.ndof]>
-                      function_space.exterior_facet_node_list)
+            self.exterior_facet_node_list = \
+                np.array(<int[:self._mesh.exterior_facets.count,:element_f.ndof]>
+                         function_space.exterior_facet_node_list)
 
         # Note that this is the function space rank and is therefore
         # always 0. The value rank may be different.
@@ -502,9 +829,9 @@
         name = "%s_nodes" % self.name
         if self._halo:
             return op2.Set(self.dof_classes, name,
-                           halo=self._halo.op2_halo)
+                           halo=self._halo.op2_halo, layers=self._mesh.layers)
         else:
-            return op2.Set(self.dof_count, name)
+            return op2.Set(self.node_count, name, layers=self._mesh.layers)
 
     @utils.cached_property
     def dof_dset(self):
@@ -526,7 +853,8 @@
         underlying mesh to the :attr:`node_set` of this :class:FunctionSpace."""
         return op2.Map(self._mesh.cell_set, self.node_set,
                        self.fiat_element.space_dimension(),
-                       self.cell_node_list - 1, "%s_cell_dof" % (self.name))
+                       self.cell_node_list - 1, "%s_cell_dof" % (self.name),
+                       offset=self.offset)
 
     @utils.cached_property
     def interior_facet_node_map(self):
@@ -558,11 +886,13 @@
     def mesh(self):
         """The :class:`Mesh` used to construct this :class:`FunctionSpace`."""
         return self._mesh
+        
+        
+        
 
 class VectorFunctionSpace(FunctionSpace):
-    def __init__(self, mesh, family, degree, dim = None, name = None):
-        super(VectorFunctionSpace, self).__init__(mesh, family, degree, name)
-
+    def __init__(self, mesh, family, degree, dim=None, name=None, vfamily=None, vdegree=None):
+        super(VectorFunctionSpace, self).__init__(mesh, family, degree, name, vfamily=vfamily, vdegree=vdegree)
         self.rank = 1
         if dim is None:
             # VectorFunctionSpace dimension defaults to the geometric dimension of the mesh.
@@ -596,7 +926,7 @@
         else:
             raise NotImplementedError("Can't make a Function defined on a "
                                       +str(type(function_space)))
-
+            
         ufl.Coefficient.__init__(self, self._function_space.ufl_element())
 
         self.name = name

=== modified file 'python/firedrake/fluidity_types.pxd'
--- python/firedrake/fluidity_types.pxd	2013-09-19 13:10:25 +0000
+++ python/firedrake/fluidity_types.pxd	2013-10-03 16:15:40 +0000
@@ -66,3 +66,4 @@
 cdef extern halo_t halo_f(void *, halo_entity)
 cdef extern void function_space_destructor_f(void *)
 cdef extern void vector_field_destructor_f(void *)
+cdef extern function_space_t extruded_mesh_f(void *, element_t *, int *)

=== modified file 'python/firedrake/io.py'
--- python/firedrake/io.py	2013-09-10 13:19:41 +0000
+++ python/firedrake/io.py	2013-10-03 16:15:40 +0000
@@ -1,6 +1,7 @@
 from evtk.hl import *
 from evtk.vtk import _get_byte_order
 from evtk.hl import _requiresLargeVTKFileSize, _addDataToFile, _appendDataToFile
+from ufl import Cell, OuterProductCell
 import numpy as np
 import os
 from pyop2.mpi import *
@@ -8,18 +9,23 @@
 
 # Dictionary used to translate the cellname of firedrake
 # to the celltype of evtk module.
-_cells = {
-    "interval": VtkLine,
-    "triangle": VtkTriangle,
-    "tetrahedron": VtkTetra
-}
-
-_points_per_cell = {
-    "interval": 2,
-    "triangle": 3,
-    "tetrahedron": 4
-}
-
+_cells = {}
+_cells[Cell("interval")] = VtkLine
+_cells[Cell("interval", 2)] = VtkLine
+_cells[Cell("interval", 3)] = VtkLine
+_cells[Cell("triangle")] = VtkTriangle
+_cells[Cell("triangle", 3)] = VtkTriangle
+_cells[Cell("tetrahedron")] = VtkTetra
+_cells[OuterProductCell(Cell("triangle"),Cell("interval"))] = VtkWedge
+
+_points_per_cell = {}
+_points_per_cell[Cell("interval")] = 2
+_points_per_cell[Cell("interval", 2)] = 2
+_points_per_cell[Cell("interval", 3)] = 2
+_points_per_cell[Cell("triangle")] = 3
+_points_per_cell[Cell("triangle", 3)] = 3
+_points_per_cell[Cell("tetrahedron")] = 4
+_points_per_cell[OuterProductCell(Cell("triangle"),Cell("interval"))] = 6
 
 class File(object):
 
@@ -114,20 +120,34 @@
         # is not 1.
 
         e = function.function_space().ufl_element()
-
         #finite element supported : CG1 and DG0
-        if not ((e.degree() == 1 and e.family() == "Lagrange") or
-                (e.degree() == 0 and e.family() == "Discontinuous Lagrange")):
+        if not ((e.family() == "Lagrange" and e.degree() == 1) or
+                (e.family() == "Discontinuous Lagrange" and e.degree() == 0) or
+                (e.family() == "OuterProductElement")):
             raise ValueError("ufl element given is not supported.")
         function.dat.halo_exchange_begin()
         function.dat.halo_exchange_end()
         mesh = function.function_space().mesh()
-        num_points = mesh._entities[0]
-        num_cells = mesh._entities[-1]
-
+        if not (e.family() == "OuterProductElement"):
+            num_points = mesh._entities[0]
+            num_cells = mesh._entities[-1]
+        else:
+            num_points = mesh._entities[0]*mesh.layers
+            num_cells = mesh._entities[-1]*(mesh.layers - 1)
+            
         # In firedrake, the indexing of the points starts from 1, instead of 0.
         # However the paraview expect them to have indexing starting from 0.
-        connectivity = mesh._cells.flatten() - 1
+        
+        if not (e.family() == "OuterProductElement"):
+            connectivity = mesh._cells.flatten() - 1
+        else:
+            # Horribly slow, hard-coded first attempt
+            connectivity_temp = np.empty([mesh.num_cells(), mesh.layers-1, _points_per_cell[mesh._ufl_cell]],dtype="int32")
+            for (ii, verts) in enumerate(mesh._cells):
+                for nl in range(mesh.layers-1):  # number of vert cells = layers - 1
+                    connectivity_temp[ii, nl] = np.concatenate((mesh.layers*(mesh._cells[ii]-1)+nl,mesh.layers*(mesh._cells[ii]-1)+nl+1))
+            connectivity = connectivity_temp.flatten() # no need to subtract 1
+ 
         data = {function.name: function.dat.data.flatten()}
 
         coordinates = self._fd_to_evtk_coord(mesh._coordinates)
@@ -135,24 +155,43 @@
         cell_types = np.empty(num_cells, dtype="uint8")
 
         # Assume that all cells are of same shape.
-        cell_types[:] = _cells[mesh._ufl_cell.cellname()].tid
-        p_c = _points_per_cell[mesh._ufl_cell.cellname()]
+        cell_types[:] = _cells[mesh._ufl_cell].tid
+        p_c = _points_per_cell[mesh._ufl_cell]
 
         # This tells which are the last nodes of each cell.
         offsets = np.arange(
             start=p_c, stop=p_c * (num_cells + 1), step=p_c, dtype='int32')
-        if e.degree() == 1:
-            large_file_flag = _requiresLargeVTKFileSize("VtkUnstructuredGrid",
-                                                        numPoints=num_points,
-                                                        numCells=num_cells,
-                                                        pointData=data,
-                                                        cellData=None)
+        if (e.family() == "Lagrange" and e.degree() == 1):     # if CG1
+            large_file_flag = _requiresLargeVTKFileSize("VtkUnstructuredGrid",
+                                                        numPoints=num_points,
+                                                        numCells=num_cells,
+                                                        pointData=data,
+                                                        cellData=None)
+        elif (e.family() == "Discontinuous Lagrange" and e.degree() == 0): # if DG0
+            large_file_flag = _requiresLargeVTKFileSize("VtkUnstructuredGrid",
+                                                        numPoints=num_points,
+                                                        numCells=num_cells,
+                                                        pointData=None,
+                                                        cellData=data)
+        elif (e.family() == "OuterProductElement" and 
+              e._A.family() == "Lagrange" and e._A.degree() == 1 and 
+              e._B.family() == "Lagrange" and e._B.degree() == 1):
+            large_file_flag = _requiresLargeVTKFileSize("VtkUnstructuredGrid",
+                                                        numPoints=num_points,
+                                                        numCells=num_cells,
+                                                        pointData=data,
+                                                        cellData=None)
+
+        elif (e.family() == "OuterProductElement" and 
+              e._A.family() == "Discontinuous Lagrange" and e._A.degree() == 0 and 
+              e._B.family() == "Discontinuous Lagrange" and e._B.degree() == 0):
+            large_file_flag = _requiresLargeVTKFileSize("VtkUnstructuredGrid",
+                                                        numPoints=num_points,
+                                                        numCells=num_cells,
+                                                        pointData=None,
+                                                        cellData=data)
         else:
-            large_file_flag = _requiresLargeVTKFileSize("VtkUnstructuredGrid",
-                                                        numPoints=num_points,
-                                                        numCells=num_cells,
-                                                        pointData=None,
-                                                        cellData=data)
+            raise ValueError("Can only support P1, P0, P1xP1, P0xP0. Hire another UROP for Summer 2014???")
 
         new_name = self._filename
 

=== modified file 'python/firedrake/mesh.py'
--- python/firedrake/mesh.py	2013-09-30 13:01:36 +0000
+++ python/firedrake/mesh.py	2013-10-03 16:15:40 +0000
@@ -37,7 +37,6 @@
 _3dexts = [".face"]
 _pexts = [".halo"]
 
-
 def _build_msh_file(input, output, dimension):
     if gmshpy:
         # We've got the gmsh python interface available, so
@@ -81,7 +80,6 @@
             else:
                 raise e
 
-
 def _get_msh_file(source, name, dimension):
 
     """Given a source code, name and dimension  of the mesh,

=== modified file 'python/firedrake/solving.py'
--- python/firedrake/solving.py	2013-10-03 14:18:31 +0000
+++ python/firedrake/solving.py	2013-10-03 16:15:40 +0000
@@ -149,6 +149,7 @@
     2-forms. The last of these may change to a native Firedrake type
     in a future release.
     """
+
     kernels = ffc_interface.compile_form(f, "form")
 
     fd = f.form_data()
@@ -163,6 +164,7 @@
 
     if is_mat:
         test, trial = fd.original_arguments
+        
         m = test.function_space().mesh()
         key = (compute_form_signature(f), test.function_space().dof_dset,
                trial.function_space().dof_dset,
@@ -189,6 +191,7 @@
         result = lambda : result_function
     else:
         m = coords.function_space().mesh()
+        
         # 0-forms are always scalar
         tensor = op2.Global(1, [0.0])
         result = lambda : tensor.data[0]
@@ -205,10 +208,20 @@
                                     flatten=True)
             else:
                 tensor_arg = tensor(op2.INC)
-            args = [kernel, m.cell_set, tensor_arg,
-                    coords.dat(op2.READ, coords.cell_node_map, flatten=True)]
+
+            itspace = m.cell_set
+            if itspace.layers > 1:
+                coords_xtr = m._coordinate_field
+                args = [kernel, itspace, tensor_arg,
+                        coords_xtr.dat(op2.READ, coords_xtr.cell_node_map,
+                                       flatten=True)]
+            else:
+                args = [kernel, itspace, tensor_arg,
+                        coords.dat(op2.READ, coords.cell_node_map, flatten=True)]
+
             for c in fd.original_coefficients:
                 args.append(c.dat(op2.READ, c.cell_node_map))
+
             op2.par_loop(*args)
         if domain_type == 'exterior_facet':
             if op2.MPI.parallel: raise \

=== added directory 'tests/firedrake_extrusion_assemble_p0'
=== added file 'tests/firedrake_extrusion_assemble_p0/Makefile'
--- tests/firedrake_extrusion_assemble_p0/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_assemble_p0/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_assemble_p0/firedrake_assemble_p0.xml'
--- tests/firedrake_extrusion_assemble_p0/firedrake_assemble_p0.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_assemble_p0/firedrake_assemble_p0.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_assemble_p0 as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] == "0.5"</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_assemble_p0/test_firedrake_assemble_p0.py'
--- tests/firedrake_extrusion_assemble_p0/test_firedrake_assemble_p0.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_assemble_p0/test_firedrake_assemble_p0.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,88 @@
+from firedrake import *
+import numpy as np
+import pyop2 as op2
+
+power = 5
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field. 
+# TODO: provide a kernel which will describe how coordinates are extruded.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel) 
+
+#import pyop2.configuration as cfg
+#cfg.configure(debug=1)
+
+def integrate_assemble_p0(family, degree):
+    fs = firedrake.FunctionSpace(mesh, family, degree, name="fs")
+    
+    f = firedrake.Function(fs)
+
+    fs1 = firedrake.FunctionSpace(mesh, family, degree, name="fs1")
+
+    f_rhs = firedrake.Function(fs1)
+    
+    populate_p0 = op2.Kernel("""
+void populate_tracer(double *x[], double *c[])
+{
+  x[0][0] = (c[1][2] + c[0][2]) / 2;
+}""", "populate_tracer")
+
+    coords = f.function_space().mesh()._coordinate_field
+
+    op2.par_loop(populate_p0, f.cell_set,
+            f.dat(op2.INC, f.cell_node_map),
+            coords.dat(op2.READ, coords.cell_node_map))
+
+    volume = op2.Kernel("""
+void comp_vol(double *rhs[], double *x[], double *y[])
+{
+  double area = x[0][0]*(x[2][1]-x[4][1]) + x[2][0]*(x[4][1]-x[0][1])
+               + x[4][0]*(x[0][1]-x[2][1]);
+  if (area < 0)
+    area = area * (-1.0);
+  rhs[0][0] += 0.5 * area * (x[1][2] - x[0][2]) * y[0][0];
+}""", "comp_vol")
+
+    op2.par_loop(volume, f.cell_set,
+             f_rhs.dat(op2.WRITE, f_rhs.cell_node_map),
+             coords.dat(op2.READ, coords.cell_node_map),
+             f.dat(op2.READ, f.cell_node_map)
+            )
+
+    g = op2.Global(1, data=0.0, name='g')
+
+    reduction = op2.Kernel("""
+void comp_reduction(double A[1], double *x[])
+{
+  A[0] += x[0][0];
+}""", "comp_reduction")
+
+    op2.par_loop(reduction, f_rhs.cell_set,
+             g(op2.INC),
+             f_rhs.dat(op2.READ, f_rhs.cell_node_map)
+            )
+
+    return str(g.data[0])
+
+def run_test():
+    family = "DG"
+    degree = 0
+
+    return [integrate_assemble_p0(family, degree)]
+
+if __name__=="__main__":
+
+    result = run_test()
+    for i, res in enumerate(result):
+        print "Result for extruded unit cube integration "+`i+1`+": "+`res`

=== added directory 'tests/firedrake_extrusion_assembly'
=== added file 'tests/firedrake_extrusion_assembly/Makefile'
--- tests/firedrake_extrusion_assembly/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_assembly/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_assembly/firedrake_assembly.xml'
--- tests/firedrake_extrusion_assembly/firedrake_assembly.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_assembly/firedrake_assembly.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,22 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="error" language="python">
+import test_firedrake_assembly as test
+error = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Error inf norm. degree 1" language="python">assert error[0] &lt; 1.0e-14</test>
+    <test name="Error inf norm. degree 2" language="python">assert error[1] &lt; 1.0e-6</test>
+    <test name="Error inf norm. degree 3" language="python">assert error[2] &lt; 1.0e-6</test>
+    <test name="Error inf norm. degree 4" language="python">assert error[3] &lt; 1.0e-6</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_assembly/test_firedrake_assembly.py'
--- tests/firedrake_extrusion_assembly/test_firedrake_assembly.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_assembly/test_firedrake_assembly.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,57 @@
+from firedrake import *
+import numpy as np
+
+power = 1
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel)
+
+#import pyop2.configuration as cfg
+#cfg.configure(debug=1)
+
+def identity_xtr(family, degree):
+    fs = firedrake.FunctionSpace(mesh, family, degree, name="fs")
+
+    f = firedrake.Function(fs)
+    out = firedrake.Function(fs)
+
+    u = firedrake.TrialFunction(fs)
+    v = firedrake.TestFunction(fs)
+
+    M = firedrake.assemble(firedrake.dot(firedrake.grad(u),firedrake.grad(v))*firedrake.dx)
+    
+    f.interpolate(firedrake.Expression("6.0"))
+
+    b = firedrake.assemble(f*v*firedrake.dx)
+
+    firedrake.solve(u*v*firedrake.dx==f*v*firedrake.dx, out)
+
+    return np.max(np.abs(out.dat.data-f.dat.data))
+
+def run_test():
+    family = "Lagrange"
+    degree = range(1,5)
+
+    error = []
+    for d in degree:
+        error.append(identity_xtr(family, d))
+    return error
+
+
+if __name__=="__main__":
+
+    error = run_test()
+    for i, err in enumerate(error):
+        print "Inf norm error in identity for Lagrange "+`i+1`+": "+`err`

=== added directory 'tests/firedrake_extrusion_helmholtz'
=== added file 'tests/firedrake_extrusion_helmholtz/Makefile'
--- tests/firedrake_extrusion_helmholtz/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_helmholtz/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_helmholtz/test_firedrake_extrusion_helmholtz.py'
--- tests/firedrake_extrusion_helmholtz/test_firedrake_extrusion_helmholtz.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_helmholtz/test_firedrake_extrusion_helmholtz.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,59 @@
+"""This demo program solves Helmholtz's equation
+ 
+  - div grad u(x, y) + u(x,y) = f(x, y)
+ 
+on the unit square with source f given by
+ 
+  f(x, y) = (1.0 + 8.0*pi**2)*cos(x[0]*2*pi)*cos(x[1]*2*pi)
+ 
+and the analytical solution
+ 
+  u(x, y) = cos(x[0]*2*pi)*cos(x[1]*2*pi)
+"""
+  
+# Begin demo
+from firedrake import *
+import sys
+
+power = int(sys.argv[1])
+# Create mesh and define function space
+m = UnitSquareMesh(2**power,2**power)
+layers = 2**int(sys.argv[1])+1
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = %(height)s*j[0][0];
+}""" % { 'height' : str(1.0/2**int(sys.argv[1])) }
+
+mesh = ExtrudedMesh(m, layers, extrusion_kernel)
+V = FunctionSpace(mesh, "Lagrange", 1, vfamily="Lagrange", vdegree=1)
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Function(V)
+f.interpolate(Expression("(1+12*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2)*cos(x[2]*pi*2)"))
+a = (dot(grad(v),grad(u))+v*u)*dx
+L = f*v*dx
+
+# Compute solution
+A = assemble(a)
+b = assemble(L)
+x = Function(V)
+solve(a==L, x)
+
+# Analytical solution
+f.interpolate( Expression("cos(x[0]*pi*2)*cos(x[1]*pi*2)*cos(x[2]*pi*2)"))
+
+print sqrt(assemble((x-f)*(x-f)*dx))
+ 
+# Save solution in VTK format
+#file = File("helmholtz.pvd")
+#file << x
+#file << f

=== added directory 'tests/firedrake_extrusion_identity'
=== added file 'tests/firedrake_extrusion_identity/Makefile'
--- tests/firedrake_extrusion_identity/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_identity/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_identity/firedrake_extrusion.xml'
--- tests/firedrake_extrusion_identity/firedrake_extrusion.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_identity/firedrake_extrusion.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,22 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="error" language="python">
+import test_firedrake_extrusion as test
+error = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Error inf norm. degree 1" language="python">assert error[0] &lt; 1.0e-14</test>
+    <test name="Error inf norm. degree 2" language="python">assert error[1] &lt; 1.0e-6</test>
+    <test name="Error inf norm. degree 3" language="python">assert error[2] &lt; 1.0e-6</test>
+    <test name="Error inf norm. degree 4" language="python">assert error[3] &lt; 1.0e-6</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_identity/test_firedrake_extrusion.py'
--- tests/firedrake_extrusion_identity/test_firedrake_extrusion.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_identity/test_firedrake_extrusion.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,58 @@
+from firedrake import *
+import numpy as np
+
+power = 1
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field.
+# TODO: provide a kernel which will describe how coordinates are extruded.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel)
+
+#import pyop2.configuration as cfg
+#cfg.configure(debug=1)
+
+def identity_xtr(family, degree):
+    fs = firedrake.FunctionSpace(mesh, family, degree, name="fs")
+
+    f = firedrake.Function(fs)
+    out = firedrake.Function(fs)
+
+    u = firedrake.TrialFunction(fs)
+    v = firedrake.TestFunction(fs)
+
+    M = firedrake.assemble(u*v*firedrake.dx)
+
+    f.interpolate(firedrake.Expression("x[0]"))
+
+    b = firedrake.assemble(f*v*firedrake.dx)
+
+    firedrake.solve(u*v*firedrake.dx==f*v*firedrake.dx, out)
+
+    return np.max(np.abs(out.dat.data-f.dat.data))
+
+def run_test():
+    family = "Lagrange"
+    degree = range(1,5)
+
+    error = []
+    for d in degree:
+        error.append(identity_xtr(family, d))
+    return error
+
+
+if __name__=="__main__":
+
+    error = run_test()
+    for i, err in enumerate(error):
+        print "Inf norm error in identity for Lagrange "+`i+1`+": "+`err`

=== added directory 'tests/firedrake_extrusion_lsw'
=== added file 'tests/firedrake_extrusion_lsw/Makefile'
--- tests/firedrake_extrusion_lsw/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_lsw/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_lsw/firedrake_extrusion_lsw.xml'
--- tests/firedrake_extrusion_lsw/firedrake_extrusion_lsw.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_lsw/firedrake_extrusion_lsw.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="special" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_extrusion_lsw as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] &lt; 0.0002</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_lsw/test_firedrake_extrusion_lsw.py'
--- tests/firedrake_extrusion_lsw/test_firedrake_extrusion_lsw.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_lsw/test_firedrake_extrusion_lsw.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,77 @@
+"""Demo of the Shallow Water
+"""
+  
+# Begin demo
+from firedrake import *
+ 
+power = 5
+# Create mesh and define function space
+m = UnitSquareMesh(2**power,2**power)
+layers = 5
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.25*j[0][0];
+}"""
+
+mesh = ExtrudedMesh(m, layers, extrusion_kernel)
+W = FunctionSpace(mesh, "BDM", 1, vfamily="Lagrange", vdegree=1)
+X = FunctionSpace(mesh, "DG", 0, vfamily="Lagrange", vdegree=1)
+Xplot = FunctionSpace(mesh, "CG", 1, vfamily="Lagrange", vdegree=1)
+
+# Define starting field
+u_0 = Function(W)
+u_h = Function(W)
+u_1 = Function(W)
+p_0 = Function(X)
+p_1 = Function(X)
+p_plot = Function(Xplot)
+p_0.interpolate(Expression("sin(4*pi*x[0])*sin(2*pi*x[1])"))
+
+T = 0.5
+t = 0
+dt = 0.0025
+
+file = File("lsw3d.pvd")
+p_trial = TrialFunction(Xplot)
+p_test = TestFunction(Xplot)
+solve(p_trial*p_test*dx == p_0*p_test*dx, p_plot)
+file << p_plot, t
+
+while t < T:
+	u = TrialFunction(W)
+	w = TestFunction(W)
+	a_1 = dot(w, u)*dx
+	L_1 = dot(w, u_0)*dx + 0.5*dt*div(w)*p_0*dx
+	solve(a_1==L_1, u_h)
+
+	p = TrialFunction(X)
+	phi = TestFunction(X)
+	a_2 = phi*p*dx
+	L_2 = phi*p_0*dx - dt*phi*div(u_h)*dx
+	solve(a_2==L_2, p_1)
+
+	u = TrialFunction(W)
+	w = TestFunction(W)
+	a_3 = dot(w, u)*dx
+	L_3 = dot(w, u_h)*dx + 0.5*dt*div(w)*p_1*dx
+	solve(a_3==L_3, u_1)
+
+	#u_0.assign(u_1)
+	#p_0.assign(p_1)
+	u_0 = u_1
+	p_0 = p_1
+	t += dt
+
+	# project into P1 x P1 for plotting
+	p_trial = TrialFunction(Xplot)
+	p_test = TestFunction(Xplot)
+	solve(p_trial*p_test*dx == p_0*p_test*dx, p_plot)
+	file << p_plot, t
+	print t
\ No newline at end of file

=== added directory 'tests/firedrake_extrusion_p0'
=== added file 'tests/firedrake_extrusion_p0/Makefile'
--- tests/firedrake_extrusion_p0/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_p0/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_p0/firedrake_p0.xml'
--- tests/firedrake_extrusion_p0/firedrake_p0.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_p0/firedrake_p0.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_p0 as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] == "3.0"</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_p0/test_firedrake_p0.py'
--- tests/firedrake_extrusion_p0/test_firedrake_p0.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_p0/test_firedrake_p0.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,72 @@
+from firedrake import *
+import numpy as np
+import pyop2 as op2
+
+power = 5
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field. 
+# TODO: provide a kernel which will describe how coordinates are extruded.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel) 
+
+#import pyop2.configuration as cfg
+#cfg.configure(debug=1)
+
+def integrate_p0(family, degree):
+    fs = firedrake.FunctionSpace(mesh, family, degree, name="fs")
+    
+    f = firedrake.Function(fs)
+    
+    populate_p0 = op2.Kernel("""
+void populate_tracer(double *x[])
+{
+  x[0][0] = 3;
+}""", "populate_tracer")
+
+    op2.par_loop(populate_p0, f.cell_set,
+            f.dat(op2.INC, f.cell_node_map))
+
+    volume = op2.Kernel("""
+void comp_vol(double A[1], double *x[], double *y[])
+{
+  double area = x[0][0]*(x[2][1]-x[4][1]) + x[2][0]*(x[4][1]-x[0][1])
+               + x[4][0]*(x[0][1]-x[2][1]);
+  if (area < 0)
+    area = area * (-1.0);
+  A[0] += 0.5 * area * (x[1][2] - x[0][2]) * y[0][0];
+}""", "comp_vol")
+
+    g = op2.Global(1, data=0.0, name='g')
+
+    coords = f.function_space().mesh()._coordinate_field
+
+    op2.par_loop(volume, f.cell_set,
+             g(op2.INC),
+             coords.dat(op2.READ, coords.cell_node_map),
+             f.dat(op2.READ, f.cell_node_map)
+            )
+
+    return str(g.data[0])
+
+def run_test():
+    family = "DG"
+    degree = 0
+
+    return [integrate_p0(family, degree)]
+
+if __name__=="__main__":
+
+    result = run_test()
+    for i, res in enumerate(result):
+        print "Result for extruded unit cube integration "+`i+1`+": "+`res`

=== added directory 'tests/firedrake_extrusion_project'
=== added file 'tests/firedrake_extrusion_project/Makefile'
--- tests/firedrake_extrusion_project/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_project/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_project/test_firedrake_extrusion_project.py'
--- tests/firedrake_extrusion_project/test_firedrake_extrusion_project.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_project/test_firedrake_extrusion_project.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,45 @@
+"""This demo program projects an analytic expression into a function space
+"""
+  
+# Begin demo
+from firedrake import *
+import sys
+
+power = int(sys.argv[1])
+# Create mesh and define function space
+m = UnitSquareMesh(2**power,2**power)
+layers = 2**int(sys.argv[1])+1
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = %(height)s*j[0][0];
+}""" % { 'height' : str(1.0/2**int(sys.argv[1])) }
+
+mesh = ExtrudedMesh(m, layers, extrusion_kernel)
+V = FunctionSpace(mesh, "Lagrange", 1, vfamily="DG", vdegree=0)
+W = FunctionSpace(mesh, "Lagrange", 3, vfamily="Lagrange", vdegree=3)
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Function(V)
+g = Function(W)
+f.interpolate(Expression("cos(x[0]*pi*2)*cos(x[1]*pi*2)*cos(x[2]*pi*2)"))
+g.interpolate(Expression("cos(x[0]*pi*2)*cos(x[1]*pi*2)*cos(x[2]*pi*2)"))
+a = u*v*dx
+L = f*v*dx
+x = Function(V)
+solve(a==L, x)
+
+print sqrt(assemble((x-g)*(x-g)*dx))
+ 
+# Save solution in VTK format
+#file = File("helmholtz.pvd")
+#file << x
+#file << f
\ No newline at end of file

=== added directory 'tests/firedrake_extrusion_rhs'
=== added file 'tests/firedrake_extrusion_rhs/Makefile'
--- tests/firedrake_extrusion_rhs/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_rhs/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_rhs/firedrake_rhs.xml'
--- tests/firedrake_extrusion_rhs/firedrake_rhs.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_rhs/firedrake_rhs.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_rhs as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] == "0.5"</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_rhs/test_firedrake_rhs.py'
--- tests/firedrake_extrusion_rhs/test_firedrake_rhs.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_rhs/test_firedrake_rhs.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,62 @@
+from firedrake import *
+import numpy as np
+import pyop2 as op2
+import ufl
+
+power = 5
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field. 
+# TODO: provide a kernel which will describe how coordinates are extruded.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel) 
+
+#import pyop2.configuration as cfg
+#cfg.configure(debug=1)
+
+def integrate_rhs(family, degree):
+    horiz = ufl.FiniteElement(family, None, degree)
+    vert = ufl.FiniteElement(family, None, degree)
+    prod = ufl.OuterProductElement(horiz, vert)
+
+    fs = firedrake.FunctionSpace(mesh, prod, name="fs")
+    
+    f = firedrake.Function(fs)
+    
+    populate_p0 = op2.Kernel("""
+void populate_tracer(double *x[], double *c[])
+{
+  x[0][0] = ((c[1][2] + c[0][2]) / 2);
+}""", "populate_tracer")
+
+    coords = f.function_space().mesh()._coordinate_field
+
+    op2.par_loop(populate_p0, f.cell_set,
+            f.dat(op2.INC, f.cell_node_map),
+            coords.dat(op2.READ, coords.cell_node_map))
+
+    g = firedrake.assemble(f*firedrake.dx)
+
+    return str(g)
+
+def run_test():
+    family = "DG"
+    degree = 0
+
+    return [integrate_rhs(family, degree)]
+
+if __name__=="__main__":
+
+    result = run_test()
+    for i, res in enumerate(result):
+        print "Result for extruded unit cube integration "+`i+1`+": "+`res`

=== added directory 'tests/firedrake_extrusion_two_step'
=== added file 'tests/firedrake_extrusion_two_step/Makefile'
--- tests/firedrake_extrusion_two_step/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_two_step/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_two_step/firedrake_extrusion_two_step.xml'
--- tests/firedrake_extrusion_two_step/firedrake_extrusion_two_step.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_two_step/firedrake_extrusion_two_step.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>firedrake</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_extrusion_two_step as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] &lt; 0.0002</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_two_step/test_firedrake_extrusion_two_step.py'
--- tests/firedrake_extrusion_two_step/test_firedrake_extrusion_two_step.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_two_step/test_firedrake_extrusion_two_step.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,84 @@
+"""Testing RT stuff
+"""
+  
+# Begin demo
+from firedrake import *
+ 
+power = 4
+# Create mesh and define function space
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = ExtrudedMesh(m, layers, extrusion_kernel)
+
+def two_step():
+	V = FunctionSpace(mesh, "Lagrange", 2, vfamily="DG", vdegree=0)
+	W = FunctionSpace(mesh, "BDM", 1, vfamily="DG", vdegree=0)
+	X = FunctionSpace(mesh, "DG", 0, vfamily="DG", vdegree=0)
+
+	# Define starting field
+	f0 = Function(V)
+	f0.interpolate(Expression("1 + x[0]*x[0] + x[1]*x[1]"))
+
+	# DO IN ONE STEP
+	u = TrialFunction(X)
+	v = TestFunction(X)
+	a = u*v*dx
+	L = div(grad(f0))*v*dx
+
+	A = assemble(a)
+	b = assemble(L)
+	f_e = Function(X)
+	solve(a==L, f_e)
+	#print "One step:"
+	#print f_e.dat.data
+
+	# DO IN TWO STEPS
+	u = TrialFunction(W)
+	v = TestFunction(W)
+	a = dot(u,v)*dx
+	L = dot(grad(f0),v)*dx
+
+	# Compute solution
+	A = assemble(a)
+	b = assemble(L)
+	f1 = Function(W)
+	solve(a==L, f1)
+	# x should be (2x, 2y) but we have no way of checking............
+
+	u = TrialFunction(X)
+	v = TestFunction(X)
+	a = u*v*dx
+	L = div(f1)*v*dx
+
+	# Compute solution
+	A = assemble(a)
+	b = assemble(L)
+	f2 = Function(X)
+	solve(a==L, f2)
+
+	#print "Two steps:"
+	#print f2.dat.data
+
+	return sum(f2.dat.data - f_e.dat.data)
+
+
+def run_test():
+    return [two_step()]
+
+if __name__=="__main__":
+
+    result = run_test()
+    for i, res in enumerate(result):
+        print "Result for extruded two step error "+`i+1`+": "+`res`
\ No newline at end of file

=== added directory 'tests/firedrake_extrusion_unit_cube'
=== added file 'tests/firedrake_extrusion_unit_cube/Makefile'
--- tests/firedrake_extrusion_unit_cube/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_unit_cube/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_unit_cube/firedrake_unit_cube.xml'
--- tests/firedrake_extrusion_unit_cube/firedrake_unit_cube.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_unit_cube/firedrake_unit_cube.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>ufl</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_unit_cube as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] == "1.0"</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_unit_cube/test_firedrake_unit_cube.py'
--- tests/firedrake_extrusion_unit_cube/test_firedrake_unit_cube.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_unit_cube/test_firedrake_unit_cube.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,64 @@
+from firedrake import *
+import numpy as np
+import pyop2 as op2
+
+power = 5
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+# coordinates as a field. 
+# A kernel which describes how coordinates are extruded.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel)
+
+def integrate_unit_cube(family, degree):
+    fs = firedrake.FunctionSpace(mesh, family, degree, name="fs")
+    
+    f = firedrake.Function(fs)
+
+    volume = op2.Kernel("""
+void comp_vol(double A[1], double *x[], double *y[])
+{
+  double area = x[0][0]*(x[2][1]-x[4][1]) + x[2][0]*(x[4][1]-x[0][1])
+               + x[4][0]*(x[0][1]-x[2][1]);
+  if (area < 0)
+    area = area * (-1.0);
+  A[0] += 0.5 * area * (x[1][2] - x[0][2]);
+}""", "comp_vol")
+
+    g = op2.Global(1, data=0.0, name='g')
+
+    coords = f.function_space().mesh()._coordinate_field
+
+    op2.par_loop(volume, f.cell_set,
+             g(op2.INC),
+             coords.dat(op2.READ, coords.cell_node_map),
+             f.dat(op2.READ, f.cell_node_map)
+            )
+
+    return str(g.data[0])
+
+def run_test():
+    family = "Lagrange"
+    degree = range(1,2)
+
+    result = []
+    for d in degree:
+        result.append(integrate_unit_cube(family, d))
+    return result
+
+
+if __name__=="__main__":
+
+    result = run_test()
+    for i, res in enumerate(result):
+        print "Result for extruded unit cube integration "+`i+1`+": "+`res`

=== added directory 'tests/firedrake_extrusion_var_p0'
=== added file 'tests/firedrake_extrusion_var_p0/Makefile'
--- tests/firedrake_extrusion_var_p0/Makefile	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_var_p0/Makefile	2013-10-03 16:15:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+
+clean:
+	rm -f  square* *.vtu *.pvd *.s *.d.1 *.stat *.log-0 *.err-0 \
+	matrixdump matrixdump.info

=== added file 'tests/firedrake_extrusion_var_p0/firedrake_var_p0.xml'
--- tests/firedrake_extrusion_var_p0/firedrake_var_p0.xml	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_var_p0/firedrake_var_p0.xml	2013-10-03 16:15:40 +0000
@@ -0,0 +1,19 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>identity</name>
+  <owner userid="gb308"/>
+  <tags>ufl</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>true</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="result" language="python">
+import test_firedrake_var_p0 as test
+result = test.run_test()
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Unit Extruded Square Integration value." language="python">assert result[0] == "0.5"</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/firedrake_extrusion_var_p0/test_firedrake_var_p0.py'
--- tests/firedrake_extrusion_var_p0/test_firedrake_var_p0.py	1970-01-01 00:00:00 +0000
+++ tests/firedrake_extrusion_var_p0/test_firedrake_var_p0.py	2013-10-03 16:15:40 +0000
@@ -0,0 +1,74 @@
+from firedrake import *
+import numpy as np
+import pyop2 as op2
+
+power = 5
+m = UnitSquareMesh(2**power,2**power)
+layers = 11
+
+# Populate the coordinates of the extruded mesh by providing the
+#coordinates as a field. 
+# TODO: provide a kernel which will describe how coordinates are extruded.
+extrusion_kernel = """
+void extrusion_kernel(double *xtr[], double *x[], int* j[])
+{
+    //Only the Z-coord is increased, the others stay the same
+    xtr[0][0] = x[0][0];
+    xtr[0][1] = x[0][1];
+    xtr[0][2] = 0.1*j[0][0];
+}"""
+
+mesh = firedrake.ExtrudedMesh(m, layers, extrusion_kernel) 
+
+#import pyop2.configuration as cfg
+#cfg.configure(debug=1)
+
+def integrate_var_p0(family, degree):
+    fs = firedrake.FunctionSpace(mesh, family, degree, name="fs")
+    
+    f = firedrake.Function(fs)
+    
+    populate_p0 = op2.Kernel("""
+void populate_tracer(double *x[], double *c[])
+{
+  //printf("%f %f \\n", c[1][2], c[0][2]);
+  x[0][0] = (c[1][2] + c[0][2]) / 2;
+}""", "populate_tracer")
+
+    coords = f.function_space().mesh()._coordinate_field
+
+    op2.par_loop(populate_p0, f.cell_set,
+            f.dat(op2.INC, f.cell_node_map),
+            coords.dat(op2.READ, coords.cell_node_map))
+
+    volume = op2.Kernel("""
+void comp_vol(double A[1], double *x[], double *y[])
+{
+  double area = x[0][0]*(x[2][1]-x[4][1]) + x[2][0]*(x[4][1]-x[0][1])
+               + x[4][0]*(x[0][1]-x[2][1]);
+  if (area < 0)
+    area = area * (-1.0);
+  A[0] += 0.5 * area * (x[1][2] - x[0][2]) * y[0][0];
+}""", "comp_vol")
+
+    g = op2.Global(1, data=0.0, name='g')
+
+    op2.par_loop(volume, f.cell_set,
+             g(op2.INC),
+             coords.dat(op2.READ, coords.cell_node_map),
+             f.dat(op2.READ, f.cell_node_map)
+            )
+
+    return str(g.data[0])
+
+def run_test():
+    family = "DG"
+    degree = 0
+
+    return [integrate_var_p0(family, degree)]
+
+if __name__=="__main__":
+
+    result = run_test()
+    for i, res in enumerate(result):
+        print "Result for extruded unit cube integration "+`i+1`+": "+`res`

=== modified file 'tests/firedrake_identity/firedrake_identity.xml'
--- tests/firedrake_identity/firedrake_identity.xml	2013-09-20 09:37:44 +0000
+++ tests/firedrake_identity/firedrake_identity.xml	2013-10-03 16:15:40 +0000
@@ -13,7 +13,7 @@
     </variable>
   </variables>
   <pass_tests>
-    <test name="Error inf norm. degree 1" language="python">assert error[0] &lt; 1.0e-14</test>
+    <test name="Error inf norm. degree 1" language="python">assert error[0] &lt; 1.0e-15</test>
     <test name="Error inf norm. degree 2" language="python">assert error[1] &lt; 1.0e-6</test>
     <test name="Error inf norm. degree 3" language="python">assert error[2] &lt; 1.0e-6</test>
     <test name="Error inf norm. degree 4" language="python">assert error[3] &lt; 1.0e-5</test>

=== modified file 'tests/firedrake_identity/test_firedrake_identity.py'
--- tests/firedrake_identity/test_firedrake_identity.py	2013-09-19 14:51:06 +0000
+++ tests/firedrake_identity/test_firedrake_identity.py	2013-10-03 16:15:40 +0000
@@ -2,7 +2,6 @@
 import numpy as np
 mesh = firedrake.UnitCubeMesh(1, 1, 1)
 
-
 def identity(family, degree):
     fs = firedrake.FunctionSpace(mesh, family, degree)
 

=== modified file 'tests/firedrake_identity_parallel/firedrake_identity_parallel.xml'
--- tests/firedrake_identity_parallel/firedrake_identity_parallel.xml	2013-09-20 09:37:44 +0000
+++ tests/firedrake_identity_parallel/firedrake_identity_parallel.xml	2013-10-03 16:15:40 +0000
@@ -14,7 +14,7 @@
     </variable>
   </variables>
   <pass_tests>
-    <test name="Error inf norm. degree 1" language="python">assert error[0] &lt; 1.0e-14</test>
+    <test name="Error inf norm. degree 1" language="python">assert error[0] &lt; 1.0e-15</test>
     <test name="Error inf norm. degree 2" language="python">assert error[1] &lt; 1.0e-6</test>
     <test name="Error inf norm. degree 3" language="python">assert error[2] &lt; 1.0e-6</test>
     <test name="Error inf norm. degree 4" language="python">assert error[3] &lt; 1.0e-5</test>