Merge lp:~btollit/fluidity/pressure_CV_weak_pressure_BC into lp:fluidity

[Brendan Tollit]

This merge will add a procedure to be able to have non zero pressure Dirichlet boundary conditions when using a CV discretisation for pressure. As for the CG pressure BCs the integrals are included for both weak and strong Dirichlet. For multiphase the phase volume fraction is included via a FE spatial interpolation, akin to their inclusion in interior domain terms.

Code Additions:

- assemble/Pressure_Dirichlet_BCS_CV.F90. Adds the BC integrals via forming the same local matrix as the weak velocity BCs but uses it transposed, which is similar to pressure CG version(s).

Code Changes:

- assemble/Momentum_Equation.F90. If cv_pressure calls the above routine after assembling the main momentum equation terms.

- assemble/Makefile.in. Update for new file/module.

- assemble/Makefile.dependencies. Update for new file/module.

Test Cases Removed:

- tests/mphase_test1_2d. Has no .xml and appears to have not been changed in a long time. Uses porous_media options and was found when looking for suitable pressure BC test cases.

Test Cases Added:

- tests/darcy_p0p1cv_pressBCinlet. Simple channel flow model using darcy flow with a pressure gradient using p0p1cv for 1,2 and 3d. Analytic and very fast.

- tests/darcy_p0p1cv_pressBCinlet_1solidphase. Same as above but with an extra phase to represent the porous media for darcy flow. Analytic and very fast.

- tests/darcy_p0p1cv_pressBCinlet_rotated_domain. Similar to the above but with a rotated domain as the pressure BCs have a normal component so worth checking non axis aligned integrals. Analytic and very fast.

- tests/darcy_p0p1_test_cty_cv_pressBCinlet. Similar to above but with CG pressure and CV continuity. This is added as not many CG pressure BC test cases appear to exist. Analytic and very fast.

- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain. Same as above but rotated. Analytic and very fast.

- tests/darcy_p0p1_test_cty_cv_velBCinlet. Similar to above but with weak velocity BC instead. Just to test a combination of two previously merged branches and the test case is very similar to above.

Buildbot green (on Nov 25 18:48)

[lp:~btollit/fluidity/pressure_CV_weak_pressure_BC]

3872. By Brendan Tollit on 2011-12-19

Very small change to a comment.

3873. By Brendan Tollit on 2011-12-19

Merge trunk into branch

[Brendan Tollit (btollit)]

Since original merge request, added a very small change to a comment.

Merged latest trunk into branch and re-ran branch buildbot which is green.

Needs a volunteer to review please!

[Cian Wilson (cwilson)]

Looks good. Sorry for the delay in reviewing.

One comment/suggestion:

In Makefile.in:
  Pressure_Dirichlet_BCS_CV.F90 -> Pressure_Dirichlet_BCS_CV.o

Clearly, as it still compiles this doesn't make much difference!

review: Approve

[lp:~btollit/fluidity/pressure_CV_weak_pressure_BC]

3874. By Brendan Tollit on 2011-12-23

Change .F90 to .o spotted by Cian in review.

3875. By Brendan Tollit on 2011-12-23

Merge trunk into branch

=== modified file 'assemble/Makefile.dependencies'
--- assemble/Makefile.dependencies	2011-10-27 23:11:25 +0000
+++ assemble/Makefile.dependencies	2011-12-23 12:05:25 +0000
@@ -579,7 +579,8 @@
    ../include/momentum_cg.mod ../include/momentum_dg.mod \
    ../include/momentum_diagnostic_fields.mod ../include/multiphase_module.mod \
    ../include/oceansurfaceforcing.mod ../include/parallel_tools.mod \
-   ../include/petsc_solve_state_module.mod ../include/profiler.mod \
+   ../include/petsc_solve_state_module.mod \
+   ../include/pressure_dirichlet_bcs_cv.mod ../include/profiler.mod \
    ../include/reduced_model_runtime.mod \
    ../include/rotated_boundary_conditions.mod ../include/slope_limiters_dg.mod \
    ../include/solvers.mod ../include/sparse_matrices_fields.mod \
@@ -645,6 +646,21 @@
    ../include/fdebug.h ../include/fields.mod ../include/fldebug.mod \
    ../include/spud.mod ../include/state_module.mod
 
+../include/pressure_dirichlet_bcs_cv.mod: Pressure_Dirichlet_BCS_CV.o
+	@true
+
+Pressure_Dirichlet_BCS_CV.o ../include/pressure_dirichlet_bcs_cv.mod: \
+   Pressure_Dirichlet_BCS_CV.F90 ../include/boundary_conditions.mod \
+   ../include/cv_face_values.mod ../include/cv_faces.mod \
+   ../include/cv_fields.mod ../include/cv_options.mod \
+   ../include/cv_shape_functions.mod ../include/cv_upwind_values.mod \
+   ../include/cvtools.mod ../include/diagnostic_fields.mod ../include/fdebug.h \
+   ../include/fetools.mod ../include/field_derivatives.mod \
+   ../include/field_options.mod ../include/fields.mod ../include/futils.mod \
+   ../include/global_parameters.mod ../include/multiphase_module.mod \
+   ../include/quadrature.mod ../include/sparsity_patterns_meshes.mod \
+   ../include/spud.mod ../include/state_module.mod
+
 ../include/pseudo_supermesh.mod: Pseudo_supermesh.o
 	@true
 
 
=== modified file 'assemble/Makefile.in'
--- assemble/Makefile.in	2011-10-21 10:05:25 +0000
+++ assemble/Makefile.in	2011-12-23 12:05:25 +0000
@@ -76,7 +76,7 @@
 	Linear_Shallow_Water.o \
 	Zoltan_global_variables.o Zoltan_integration.o Zoltan_callbacks.o Zoltan_detectors.o \
         Diagnostic_Children.o Reduced_Model_Runtime.o Turbine.o Implicit_Solids.o \
-				Manifold_Projections.o Hybridized_Helmholtz.o Burgers_Assembly.o
+				Manifold_Projections.o Hybridized_Helmholtz.o Burgers_Assembly.o Pressure_Dirichlet_BCS_CV.o
 
 .SUFFIXES: .F90 .c .cpp .o .a
 
 
=== modified file 'assemble/Momentum_Equation.F90'
--- assemble/Momentum_Equation.F90	2011-11-23 17:00:55 +0000
+++ assemble/Momentum_Equation.F90	2011-12-23 12:05:25 +0000
@@ -74,6 +74,7 @@
       use slope_limiters_dg
       use implicit_solids
       use multiphase_module
+      use pressure_dirichlet_bcs_cv
 
       implicit none
 
@@ -589,6 +590,11 @@
                      include_pressure_and_continuity_bcs=.not. cv_pressure)
             end if
             
+            ! If CV pressure then add in any dirichlet pressure BC integrals to the mom_rhs.
+            if (cv_pressure) then
+               call add_pressure_dirichlet_bcs_cv(mom_rhs(istate), u, p, state(istate))
+            end if
+            
             ! Add in multiphase interactions (e.g. fluid-particle drag) if necessary
             ! Note: this is done outside of construct_momentum_cg/dg to keep things
             ! neater in Momentum_CG/DG.F90, since we would need to pass around multiple phases 
 
=== added file 'assemble/Pressure_Dirichlet_BCS_CV.F90'
--- assemble/Pressure_Dirichlet_BCS_CV.F90	1970-01-01 00:00:00 +0000
+++ assemble/Pressure_Dirichlet_BCS_CV.F90	2011-12-23 12:05:25 +0000
@@ -0,0 +1,280 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!    
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    [e-mail]
+!    
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied warranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+#include "fdebug.h"
+
+module pressure_dirichlet_bcs_cv
+
+  use quadrature
+  use fields
+  use field_derivatives
+  use state_module
+  use futils
+  use fetools
+  use spud
+  use cv_faces
+  use cv_shape_functions
+  use cvtools
+  use cv_fields
+  use cv_upwind_values
+  use cv_face_values
+  use cv_options
+  use diagnostic_fields, only: calculate_diagnostic_variable
+  use boundary_conditions
+  use global_parameters, only: OPTION_PATH_LEN
+  use field_options, only: get_coordinate_field
+  use sparsity_patterns_meshes
+  use multiphase_module
+
+  implicit none
+
+  private
+  
+  public :: add_pressure_dirichlet_bcs_cv
+
+contains
+
+    ! --------------------------------------------------------------------------------------
+
+    subroutine add_pressure_dirichlet_bcs_cv(mom_rhs, u, p, state)
+       
+      !!< Add any CV pressure dirichlet BC integrals to the mom_rhs field if required.
+      !!< If this is a multiphase simulation then the phase volume fraction spatial 
+      !!< interpolation uses finite elements. This isnt striclty correct if the 
+      !!< phase volume fractions were solved with a control volume discretisation 
+      !!< and also whether of not they themselves have a weak BC is ignored. If the 
+      !!< latter was to be considered then the upwind value should be used.
+      
+      ! Note the pressure BC values come from the end of time step Pressure field
+      ! rather than using a theta average. The latter was tried but didnt work 
+      ! because OldPressure had no BC info. Taking the end of time step value 
+      ! for the BC is also what the pressure CG routine(s) do.
+      
+      type(vector_field), intent(inout) :: mom_rhs
+      type(vector_field), intent(in) :: u
+      type(scalar_field), intent(in) :: p
+      type(state_type), intent(inout) :: state ! required for phase volume fraction
+
+      ! local variables
+
+      ! degree of quadrature over cv faces
+      integer :: quaddegree
+
+      real, dimension(:,:), allocatable :: x_ele, x_ele_bdy
+      real, dimension(:,:), allocatable :: normal_bdy
+      real, dimension(:), allocatable :: detwei_bdy
+      integer, dimension(:), allocatable :: u_nodes_bdy
+
+      integer :: ele, sele, iloc, jloc, face, gi, ggi, dim
+
+      ! information about cv faces
+      type(cv_faces_type) :: cvfaces
+      ! shape functions for surface
+      type(element_type) :: x_cvbdyshape
+      type(element_type) :: u_cvbdyshape
+      
+      ! pointer to coordinates
+      type(vector_field), pointer :: x
+      
+      ! pressure BC info
+      integer, dimension(:), allocatable :: pressure_bc_type
+      real, dimension(:), allocatable :: pressure_bc_val
+      type(scalar_field) :: pressure_bc
+      
+      ! local ct matrix
+      real, dimension(:,:,:), allocatable :: ct_mat_local_bdy
+
+      ! Multiphase variables
+      logical :: multiphase
+      ! Volume fraction fields
+      type(scalar_field), pointer :: vfrac
+      type(scalar_field) :: nvfrac
+      
+      ! Variables used for FE interpolation of vfrac on surface
+      ! Values of nvfrac at the faces 
+      real, dimension(:), allocatable :: nvfrac_gi_f
+      ! CV shape functions for nvfrac (computed only if the Coordinate and
+      ! PhaseVolumeFraction meshes are different, otherwise it is
+      ! assigned x_cvbdyshape)
+      type(element_type) :: nvfrac_cvbdyshape
+
+      ewrite(1,*) 'In add_cv_pressure_weak_dirichlet_bcs'
+
+      x => extract_vector_field(state, "Coordinate")
+      
+      allocate(x_ele(x%dim, x%mesh%shape%loc))
+
+      call get_option("/geometry/quadrature/controlvolume_surface_degree", quaddegree, default=1)
+
+      cvfaces = find_cv_faces(vertices   = ele_vertices(p, 1), &
+                              dimension  = mesh_dim(p), &
+                              polydegree = p%mesh%shape%degree, &
+                              quaddegree = quaddegree)
+
+      x_cvbdyshape = make_cvbdy_element_shape(cvfaces, x%mesh%faces%shape)
+      u_cvbdyshape = make_cvbdy_element_shape(cvfaces, u%mesh%faces%shape)
+   
+      assert(surface_element_count(p) == surface_element_count(u))
+      
+      ! get the pressure BC field
+      allocate(pressure_bc_type(surface_element_count(p)))
+      
+      call get_entire_boundary_condition(p, (/"weakdirichlet", &
+                                              "dirichlet    "/), pressure_bc, pressure_bc_type)
+      
+      allocate(x_ele_bdy(x%dim,x%mesh%faces%shape%loc), &
+               detwei_bdy(x_cvbdyshape%ngi), &
+               normal_bdy(x%dim, x_cvbdyshape%ngi), &
+               pressure_bc_val(pressure_bc%mesh%shape%loc))
+      
+      allocate(u_nodes_bdy(u%mesh%faces%shape%loc))
+            
+      allocate(ct_mat_local_bdy(x%dim, p%mesh%faces%shape%loc, u%mesh%faces%shape%loc))
+
+      ! Check if we need to multiply through by the non-linear volume fraction
+      if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+         multiphase = .true.
+
+         vfrac => extract_scalar_field(state, "PhaseVolumeFraction")
+         call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction")
+         call zero(nvfrac)
+         call get_nonlinear_volume_fraction(state, nvfrac)
+         
+         assert(surface_element_count(p) == surface_element_count(vfrac))
+         
+         ewrite_minmax(nvfrac)
+         
+         ! If the Coordinate and PhaseVolumeFraction meshes are different, then we need to
+         ! generate the PhaseVolumeFraction CV shape functions.
+         if(.not.(nvfrac%mesh == x%mesh)) then
+            nvfrac_cvbdyshape = make_cvbdy_element_shape(cvfaces, nvfrac%mesh%faces%shape)
+         else
+            nvfrac_cvbdyshape = x_cvbdyshape
+            ! incease reference for deallocates below
+            call incref(nvfrac_cvbdyshape)
+         end if
+
+         allocate(nvfrac_gi_f(nvfrac_cvbdyshape%ngi))
+                  
+      else
+         multiphase = .false.
+         nullify(vfrac)
+      end if
+  
+      surface_element_loop: do sele = 1, surface_element_count(p)
+  
+        ! cycle if this not a dirichlet pressure BC.
+        if (pressure_bc_type(sele) == 0) cycle
+            
+        ele         = face_ele(x, sele)
+        x_ele       = ele_val(x, ele)
+        x_ele_bdy   = face_val(x, sele)
+        u_nodes_bdy = face_global_nodes(u, sele)
+        
+        ! Get the phase volume fraction face value 
+        if(multiphase) then
+           nvfrac_gi_f = face_val_at_quad(nvfrac, sele, nvfrac_cvbdyshape)
+        end if
+  
+        call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, &
+                              x_cvbdyshape, normal_bdy, detwei_bdy)
+  
+        ct_mat_local_bdy = 0.0
+        
+        ! calculate the ct local matrix
+        surface_nodal_loop_i: do iloc = 1, p%mesh%faces%shape%loc
+  
+           surface_face_loop: do face = 1, cvfaces%sfaces
+          
+              if(cvfaces%sneiloc(iloc,face) /= 0) then
+  
+                 surface_quadrature_loop: do gi = 1, cvfaces%shape%ngi
+  
+                    ggi = (face-1)*cvfaces%shape%ngi + gi
+                                        
+                    surface_nodal_loop_j: do jloc = 1, u%mesh%faces%shape%loc
+                       
+                       surface_inner_dimension_loop: do dim = 1, size(normal_bdy,1)
+
+                          if(multiphase) then
+                             ct_mat_local_bdy(dim, iloc, jloc) =  ct_mat_local_bdy(dim, iloc, jloc) + &
+                                   u_cvbdyshape%n(jloc,ggi)*detwei_bdy(ggi)*nvfrac_gi_f(ggi)*normal_bdy(dim, ggi)
+                          else
+                             ct_mat_local_bdy(dim, iloc, jloc) =  ct_mat_local_bdy(dim, iloc, jloc) + &
+                                   u_cvbdyshape%n(jloc,ggi)*detwei_bdy(ggi)*normal_bdy(dim, ggi)
+                          end if
+    
+                       end do surface_inner_dimension_loop
+  
+                    end do surface_nodal_loop_j
+  
+                 end do surface_quadrature_loop
+  
+              end if
+  
+           end do surface_face_loop
+  
+        end do surface_nodal_loop_i
+        
+        ! pressure dirichlet BC is -c*press_bc_val = -press_bc_val*ct, integrated over surface elements appropriate
+        surface_outer_dimension_loop: do dim = 1, size(normal_bdy,1)
+
+          ! for weak and strong pressure dirichlet bcs:
+          !      /
+          ! add -|  N_i M_j \vec n p_j, where p_j are the prescribed bc values
+          !      /
+          
+          call addto(mom_rhs, dim, u_nodes_bdy, -matmul(ele_val(pressure_bc, sele), ct_mat_local_bdy(dim,:,:)))
+  
+        end do surface_outer_dimension_loop
+  
+      end do surface_element_loop
+  
+      call deallocate(pressure_bc)
+      deallocate(pressure_bc_type)
+      deallocate(pressure_bc_val)
+      call deallocate(x_cvbdyshape)
+      call deallocate(u_cvbdyshape)
+      deallocate(x_ele_bdy, detwei_bdy, normal_bdy)
+      
+      deallocate(u_nodes_bdy)
+      deallocate(ct_mat_local_bdy)
+        
+      call deallocate(cvfaces)
+      deallocate(x_ele)
+
+      if(multiphase) then
+         call deallocate(nvfrac)
+         deallocate(nvfrac_gi_f)
+         call deallocate(nvfrac_cvbdyshape)         
+      end if
+       
+    end subroutine add_pressure_dirichlet_bcs_cv
+    
+    ! --------------------------------------------------------------------------------------
+
+end module pressure_dirichlet_bcs_cv
+
 
=== added directory 'tests/darcy_p0p1_test_cty_cv_pressBCinlet'
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile	2011-12-23 12:05:25 +0000
@@ -0,0 +1,32 @@
+SHELL = sh
+
+SIM = darcy_p0p1_test_cty_cv_pressBCinlet
+
+# options for 1d mesh interval script
+MESH_SIZE = 100.0
+LEFT_X    = 0.0
+RIGHT_X   = 300.0
+
+input: clean
+	interval --dx=$(MESH_SIZE) $(LEFT_X) $(RIGHT_X) $(SIM)_1d
+	gmsh -2 square.geo -o square.msh
+	gmsh -3 cube.geo -o cube.msh
+
+clean:
+	rm -f *.ele
+	rm -f *.node
+	rm -f *.bound
+	rm -f *.edge
+	rm -f *.face
+	rm -f *.vtu 
+	rm -f *.pvtu 
+	rm -f *.s 
+	rm -f *.stat 
+	rm -f *.log-0 
+	rm -f *.err-0 
+	rm -f *.msh
+	rm -f *.halo
+	rm -f fluidity.err*
+	rm -f fluidity.log*
+	rm -f matrixdump*
+	rm -f first_timestep_adapted_mesh*
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo	2011-12-23 12:05:25 +0000
@@ -0,0 +1,39 @@
+// define a layer variable
+lay = 5;
+
+// define a len variable
+len = 300.0;
+
+Point(1) = {0.0, 0.0, 0.0, 1.0};
+
+Extrude {len, 0.0, 0.0} {
+  Point{1}; Layers{lay}; 
+}
+
+Extrude {0.0, len, 0.0} {
+  Line{1}; Layers{lay}; 
+}
+
+Extrude {0.0, 0.0, len} {
+  Surface{5}; Layers{lay}; 
+}
+
+// bottom
+Physical Surface(1) = {5};
+
+// top
+Physical Surface(2) = {27};
+
+// left
+Physical Surface(3) = {26};
+
+// right
+Physical Surface(4) = {18};
+
+// front
+Physical Surface(5) = {14};
+
+// back
+Physical Surface(6) = {22};
+
+Physical Volume(1) = {1};
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,118 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+  <name>darcy_p0p1_test_cty_cv_pressBCinlet</name>
+  <owner userid="btollit"/>
+  <tags>flml darcy</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>
+../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml 
+../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml 
+../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml 
+    </command_line>
+    <!-- One/two/three dimensional problem for darcy flow with one region and one material using p0p1 element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flows.-->
+  </problem_definition>
+  <variables>
+    <variable name="pressure_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_1d_1.vtu")
+pressure_1d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_1d_1.vtu")
+inter_vel_1d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="max_div_vel_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_1d_1.vtu")
+max_div_vel_1d = max(v.GetScalarRange("ControlVolumeDivergence"))
+    </variable>
+    <variable name="pressure_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_2d_1.vtu")
+pressure_2d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_2d_1.vtu")
+inter_vel_2d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="max_div_vel_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_2d_1.vtu")
+max_div_vel_2d = max(v.GetScalarRange("ControlVolumeDivergence"))
+    </variable>
+    <variable name="pressure_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_3d_1.vtu")
+pressure_3d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_3d_1.vtu")
+inter_vel_3d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="max_div_vel_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_3d_1.vtu")
+max_div_vel_3d = max(v.GetScalarRange("ControlVolumeDivergence"))
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="change_P for 1d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_1d) - min(pressure_1d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 1d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_1d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_1d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="max_div_vel_1d under tolerance 1.0e-09" language="python">
+assert max_div_vel_1d &lt; 1.0e-09
+    </test>
+    <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_2d) - min(pressure_2d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 2d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_2d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="max_div_vel_2d under tolerance 1.0e-09" language="python">
+assert max_div_vel_2d &lt; 1.0e-09
+    </test>
+    <test name="change_P for 3d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_3d) - min(pressure_3d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 3d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_3d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_3d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="max_div_vel_3d under tolerance 1.0e-09" language="python">
+assert max_div_vel_3d &lt; 1.0e-09
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,411 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_1d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1 for velocity-pressure
+
+with the continuity equation tested with the control volume dual mesh
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with pressure boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 1d and one time step is done (as nothing depends on time)
+
+the left and right boundary conditions of pressure are applied strongly.
+
+For 2d and 3d a second approach is also done via solving a pressure poisson equation then calculating the 
+velocity diagnostically from this. This uses generic scalar and vector fields and python diagnostics. 
+The latter didnt work in 1d. The same vel-press element pair is used as the mixed projection formulation.</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="darcy_p0p1_test_cty_cv_pressBCinlet_1d">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="ControlVolumeDivergence" rank="0">
+      <diagnostic field_name="Velocity">
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,427 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_2d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1 for velocity-pressure
+
+with the continuity equation tested with the control volume dual mesh
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with pressure boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 2d (although this is a 1d problem) and one time step is done (as nothing depends on time)
+
+the left and right boundary conditions of pressure are applied strongly.
+
+For 2d and 3d a second approach is also done via solving a pressure poisson equation then calculating the 
+velocity diagnostically from this. This uses generic scalar and vector fields and python diagnostics. 
+The latter didnt work in 1d. The same vel-press element pair is used as the mixed projection formulation.</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="square">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">8</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="no_outflow_top_bottom">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 9</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="ControlVolumeDivergence" rank="0">
+      <diagnostic field_name="Velocity">
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,442 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_3d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1 for velocity-pressure
+
+with the continuity equation tested with the control volume dual mesh
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with pressure boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 3d (although this is a 1d problem) and one time step is done (as nothing depends on time)
+
+the left and right boundary conditions of pressure are applied strongly.
+
+For 2d and 3d a second approach is also done via solving a pressure poisson equation then calculating the 
+velocity diagnostically from this. This uses generic scalar and vector fields and python diagnostics. 
+The latter didnt work in 1d. The same vel-press element pair is used as the mixed projection formulation.</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="cube">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">4</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">3</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="no_outflow_top_bottom">
+          <surface_ids>
+            <integer_value shape="2" rank="1">1 2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <z_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </z_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="no_outflow_front_back">
+          <surface_ids>
+            <integer_value shape="2" rank="1">5 6</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="ControlVolumeDivergence" rank="0">
+      <diagnostic field_name="Velocity">
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo	2011-12-23 12:05:25 +0000
@@ -0,0 +1,22 @@
+// define a layer variable
+lay = 5;
+
+// define a len variable
+len = 300.0;
+
+Point(1) = {0.0, 0.0, 0.0, 1.0};
+
+Extrude {len, 0.0, 0.0} {
+  Point{1}; Layers{lay}; 
+}
+
+Extrude {0.0, len, 0.0} {
+  Line{1}; Layers{lay}; 
+}
+
+Physical Line(7) = {1};
+Physical Line(8) = {4};
+Physical Line(9) = {2};
+Physical Line(10) = {3};
+
+Physical Surface(11) = {5};
 
=== added directory 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain'
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile	2011-12-23 12:05:25 +0000
@@ -0,0 +1,23 @@
+SHELL = sh
+
+input: clean
+	gmsh -2 square_rotated.geo -o square_rotated.msh
+
+clean:
+	rm -f *.ele
+	rm -f *.node
+	rm -f *.bound
+	rm -f *.edge
+	rm -f *.face
+	rm -f *.vtu 
+	rm -f *.pvtu 
+	rm -f *.s 
+	rm -f *.stat 
+	rm -f *.log-0 
+	rm -f *.err-0 
+	rm -f *.msh
+	rm -f *.halo
+	rm -f fluidity.err*
+	rm -f fluidity.log*
+	rm -f matrixdump*
+	rm -f first_timestep_adapted_mesh*
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,46 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+  <name>darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain</name>
+  <owner userid="btollit"/>
+  <tags>flml darcy</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>
+../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml 
+    </command_line>
+    <!-- Two dimensional problem (actually 1d done in 2d) for darcy flow with one region and one material using p0p1_test_cty_cv element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flow. This will use a rotated gmsh domain to test applying the weak pressure BC at a non cartesian aligned axis.-->
+  </problem_definition>
+  <variables>
+    <variable name="pressure_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d_1.vtu")
+pressure_2d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d_1.vtu")
+inter_vel_2d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_2d) - min(pressure_2d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 2d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python">
+import math
+# the analytic velocity magnitude is 10.0, with an incline of 45 degrees (= pi/4.0 radians)
+analytic_vel = 10.0 * math.sin(math.pi/4.0)
+print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_2d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,372 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="square_rotated">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">5</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">8</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="no_outflow_top_bottom">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 9</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo'
--- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo	2011-12-23 12:05:25 +0000
@@ -0,0 +1,38 @@
+// define a characteristic length
+lc = 60.0;
+
+len = 300.0;
+
+// define points
+Point(1) = {0.0,0.0,0.0,lc};
+Point(2) = {len,0.0,0.0,lc};
+Point(3) = {len,len,0.0,lc};
+Point(4) = {0.0,len,0.0,lc};
+
+// define line's
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+// define line loop
+Line Loop(1) = {1,2,3,4};
+
+// define surface
+Plane Surface(1) = {1};
+
+// Rotate the domain about the z axis by 45 degrees (Pi/4)
+Rotate {{0, 0, 1}, {0, 0, 0}, Pi/4} {
+  Surface{1};
+}
+
+// Tag the lines and surface
+Physical Line(7) = {1};
+Physical Line(8) = {2};
+Physical Line(9) = {3};
+Physical Line(10) = {4};
+
+Physical Surface(11) = {1};
+
+// Make the grid regular
+Transfinite Surface"*";
 
=== added directory 'tests/darcy_p0p1_test_cty_cv_velBCinlet'
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile	2011-12-23 12:05:25 +0000
@@ -0,0 +1,32 @@
+SHELL = sh
+
+SIM = darcy_p0p1_test_cty_cv_velBCinlet
+
+# options for 1d mesh interval script
+MESH_SIZE = 100.0
+LEFT_X    = 0.0
+RIGHT_X   = 300.0
+
+input: clean
+	interval --dx=$(MESH_SIZE) $(LEFT_X) $(RIGHT_X) $(SIM)_1d
+	gmsh -2 square.geo -o square.msh
+	gmsh -3 cube.geo -o cube.msh
+
+clean:
+	rm -f *.ele
+	rm -f *.node
+	rm -f *.bound
+	rm -f *.edge
+	rm -f *.face
+	rm -f *.vtu 
+	rm -f *.pvtu 
+	rm -f *.s 
+	rm -f *.stat 
+	rm -f *.log-0 
+	rm -f *.err-0 
+	rm -f *.msh
+	rm -f *.halo
+	rm -f fluidity.err*
+	rm -f fluidity.log*
+	rm -f matrixdump*
+	rm -f first_timestep_adapted_mesh*
 
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo	2011-12-23 12:05:25 +0000
@@ -0,0 +1,39 @@
+// define a layer variable
+lay = 5;
+
+// define a len variable
+len = 300.0;
+
+Point(1) = {0.0, 0.0, 0.0, 1.0};
+
+Extrude {len, 0.0, 0.0} {
+  Point{1}; Layers{lay}; 
+}
+
+Extrude {0.0, len, 0.0} {
+  Line{1}; Layers{lay}; 
+}
+
+Extrude {0.0, 0.0, len} {
+  Surface{5}; Layers{lay}; 
+}
+
+// bottom
+Physical Surface(1) = {5};
+
+// top
+Physical Surface(2) = {27};
+
+// left
+Physical Surface(3) = {26};
+
+// right
+Physical Surface(4) = {18};
+
+// front
+Physical Surface(5) = {14};
+
+// back
+Physical Surface(6) = {22};
+
+Physical Volume(1) = {1};
 
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,94 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+  <name>darcy_p0p1_test_cty_cv_velBCinlet</name>
+  <owner userid="btollit"/>
+  <tags>flml darcy</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>
+../../bin/fluidity darcy_p0p1_test_cty_cv_velBCinlet_1d.flml 
+../../bin/fluidity darcy_p0p1_test_cty_cv_velBCinlet_2d.flml 
+../../bin/fluidity darcy_p0p1_test_cty_cv_velBCinlet_3d.flml 
+    </command_line>
+    <!-- One/two/three dimensional problem for darcy flow with one region and one material using p0p1_test_cty_cv element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flow-->
+  </problem_definition>
+  <variables>
+    <variable name="pressure_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_1d_1.vtu")
+pressure_1d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_1d_1.vtu")
+inter_vel_1d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="pressure_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_2d_1.vtu")
+pressure_2d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_2d_1.vtu")
+inter_vel_2d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="pressure_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_3d_1.vtu")
+pressure_3d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_3d_1.vtu")
+inter_vel_3d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="change_P for 1d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_1d) - min(pressure_1d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 1d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_1d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_1d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_2d) - min(pressure_2d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 2d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_2d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="change_P for 3d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_3d) - min(pressure_3d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 3d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_3d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_3d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,391 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_velBCinlet_1d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1_test_cty_cv for velocity-pressure
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with velocity boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 1d and one time step is done (as nothing depends on time)</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="darcy_p0p1_test_cty_cv_velBCinlet_1d">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">2</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">5.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,407 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_velBCinlet_2d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1_test_cty_cv for velocity-pressure
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with velocity boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 2d (although this is a 1d problem) and one time step is done (as nothing depends on time)</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="square">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">2</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">8</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">5.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="no_outflow_top_bottom">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 9</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,404 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1_test_cty_cv_velBCinlet_3d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1_test_cty_cv for velocity-pressure
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with velocity boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 3d (although this is a 1d problem) and one time step is done (as nothing depends on time)</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="cube">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">2</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <test_continuity_with_cv_dual/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">4</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">3</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">5.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="no_outflow_top_bottom">
+          <surface_ids>
+            <integer_value shape="2" rank="1">1 2</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <boundary_conditions name="no_outflow_front_back">
+          <surface_ids>
+            <integer_value shape="2" rank="1">5 6</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo'
--- tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo	2011-12-23 12:05:25 +0000
@@ -0,0 +1,22 @@
+// define a layer variable
+lay = 5;
+
+// define a len variable
+len = 300.0;
+
+Point(1) = {0.0, 0.0, 0.0, 1.0};
+
+Extrude {len, 0.0, 0.0} {
+  Point{1}; Layers{lay}; 
+}
+
+Extrude {0.0, len, 0.0} {
+  Line{1}; Layers{lay}; 
+}
+
+Physical Line(7) = {1};
+Physical Line(8) = {4};
+Physical Line(9) = {2};
+Physical Line(10) = {3};
+
+Physical Surface(11) = {5};
 
=== added directory 'tests/darcy_p0p1cv_pressBCinlet'
=== added file 'tests/darcy_p0p1cv_pressBCinlet/Makefile'
--- tests/darcy_p0p1cv_pressBCinlet/Makefile	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1cv_pressBCinlet/Makefile	2011-12-23 12:05:25 +0000
@@ -0,0 +1,32 @@
+SHELL = sh
+
+SIM = darcy_p0p1cv_pressBCinlet
+
+# options for 1d mesh interval script
+MESH_SIZE = 100.0
+LEFT_X    = 0.0
+RIGHT_X   = 300.0
+
+input: clean
+	interval --dx=$(MESH_SIZE) $(LEFT_X) $(RIGHT_X) $(SIM)_1d
+	gmsh -2 square.geo -o square.msh
+	gmsh -3 cube.geo -o cube.msh
+
+clean:
+	rm -f *.ele
+	rm -f *.node
+	rm -f *.bound
+	rm -f *.edge
+	rm -f *.face
+	rm -f *.vtu 
+	rm -f *.pvtu 
+	rm -f *.s 
+	rm -f *.stat 
+	rm -f *.log-0 
+	rm -f *.err-0 
+	rm -f *.msh
+	rm -f *.halo
+	rm -f fluidity.err*
+	rm -f fluidity.log*
+	rm -f matrixdump*
+	rm -f first_timestep_adapted_mesh*
 
=== added file 'tests/darcy_p0p1cv_pressBCinlet/cube.geo'
--- tests/darcy_p0p1cv_pressBCinlet/cube.geo	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1cv_pressBCinlet/cube.geo	2011-12-23 12:05:25 +0000
@@ -0,0 +1,39 @@
+// define a layer variable
+lay = 5;
+
+// define a len variable
+len = 300.0;
+
+Point(1) = {0.0, 0.0, 0.0, 1.0};
+
+Extrude {len, 0.0, 0.0} {
+  Point{1}; Layers{lay}; 
+}
+
+Extrude {0.0, len, 0.0} {
+  Line{1}; Layers{lay}; 
+}
+
+Extrude {0.0, 0.0, len} {
+  Surface{5}; Layers{lay}; 
+}
+
+// bottom
+Physical Surface(1) = {5};
+
+// top
+Physical Surface(2) = {27};
+
+// left
+Physical Surface(3) = {26};
+
+// right
+Physical Surface(4) = {18};
+
+// front
+Physical Surface(5) = {14};
+
+// back
+Physical Surface(6) = {22};
+
+Physical Volume(1) = {1};
 
=== added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml'
--- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,94 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+  <name>darcy_p0p1cv_pressBCinlet</name>
+  <owner userid="btollit"/>
+  <tags>flml darcy</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>
+../../bin/fluidity darcy_p0p1cv_pressBCinlet_1d.flml 
+../../bin/fluidity darcy_p0p1cv_pressBCinlet_2d.flml 
+../../bin/fluidity darcy_p0p1cv_pressBCinlet_3d.flml 
+    </command_line>
+    <!-- One/two/three dimensional problem for darcy flow with one region and one material using p0p1cv element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flow. This also tests the use of pressure BC for CV pressure.-->
+  </problem_definition>
+  <variables>
+    <variable name="pressure_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_1d_1.vtu")
+pressure_1d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_1d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_1d_1.vtu")
+inter_vel_1d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="pressure_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_2d_1.vtu")
+pressure_2d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_2d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_2d_1.vtu")
+inter_vel_2d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+    <variable name="pressure_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_3d_1.vtu")
+pressure_3d = v.GetScalarRange("Pressure")
+    </variable>
+    <variable name="inter_vel_3d" language="python">
+import vtktools
+v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_3d_1.vtu")
+inter_vel_3d = max(v.GetScalarRange("interstitial_velocity"))
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="change_P for 1d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_1d) - min(pressure_1d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 1d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_1d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_1d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_2d) - min(pressure_2d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 2d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_2d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+    <test name="change_P for 3d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python">
+change_P = abs(max(pressure_3d) - min(pressure_3d))
+visc = 1.0e-04
+darcy_vel_BC = 5.0
+perm = 1.0e-10
+domain_length = 300.0
+print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm))
+assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) &lt; 1.0e-09
+    </test>
+    <test name="interstitial velocity 3d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python">
+analytic_vel = 10.0
+print 'Relative error of interstitial velocity: ',abs((inter_vel_3d - analytic_vel)/analytic_vel)
+assert abs((inter_vel_3d - analytic_vel)/analytic_vel) &lt; 1.0e-09
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
 
=== added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml'
--- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,384 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1cv_pressBCinlet_1d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1cv for velocity-pressure
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with pressure boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 1d and one time step is done (as nothing depends on time).</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="darcy_p0p1cv_pressBCinlet_1d">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <control_volumes/>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml'
--- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,392 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1cv_pressBCinlet_2d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1cv for velocity-pressure
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with pressure boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 2d (although this is a 1d problem) and one time step is done (as nothing depends on time).</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="square">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <control_volumes/>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">8</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="no_outflow_top_bottom">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 9</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <vector_field name="Absorption" rank="1">
+          <diagnostic>
+            <mesh name="DGP0"/>
+            <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+              <string_value lines="20" type="python"># get the prescribed fields 
+perm = states["fluid"].scalar_fields["Permeability"]
+visc = states["fluid"].scalar_fields["Viscosity"]
+
+# loop the DMmesh_p0 nodes (which are element centred)
+for n in range(field.node_count):
+   perm_n = perm.node_val(n)
+   visc_n = visc.node_val(n)
+   
+   # calc the absorption term
+   sigma_n = visc_n/perm_n
+   
+   # set the absorption term
+   field.set(n,sigma_n)</string_value>
+              <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+            </algorithm>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+          <include_pressure_correction/>
+        </vector_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="Porosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.5</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Permeability" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-10</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Viscosity" rank="0">
+      <prescribed>
+        <mesh name="DGP0"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-04</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="interstitial_velocity" rank="1">
+      <diagnostic>
+        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
+          <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"]
+darcy_vel = states["fluid"].vector_fields["Velocity"]
+
+for ele in range(field.element_count):
+   phi_ele = phi.node_val(ele)
+   
+   factor_ele = 1.0/max(phi_ele,1.0e-08)
+   
+   vel_ele = []
+   for i in range(len(darcy_vel.ele_val(ele))):
+     vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i])
+   
+   field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value>
+          <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements
+
+also assume that the velocity is dg</comment>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml'
--- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml	1970-01-01 00:00:00 +0000
+++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml	2011-12-23 12:05:25 +0000
@@ -0,0 +1,398 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">darcy_p0p1cv_pressBCinlet_3d</string_value>
+    <comment>a simple short test case for darcy flow using the element type p0p1cv for velocity-pressure
+
+the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution)
+
+it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct 
+
+the darcy flow equation is
+
+sigma*darcy_vel = - grad P
+
+sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options)
+
+python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity
+
+this case has a 1 region 1 material with pressure boundary condition inflow. 
+
+to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input.
+
+the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) 
+
+the geometry is 3d (although this is a 1d problem) and one time step is done (as nothing depends on time).</comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="cube">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DGP0">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="OutputMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period_in_timesteps>
+      <constant>
+        <integer_value rank="0">0</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="OutputMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <control_volumes/>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-10</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="right_outflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">4</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="left_inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">3</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <constant>
+              <real_value rank="0">1.5e+09</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <none/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>