Merge lp:~fluidity-core/fluidity/gmsh-on-sphere into lp:fluidity

[Stephan Kramer]

This merge fixes the reading of gmsh input meshes in fluidity, where the input mesh is 2d spherical, to be extruded inside fluidity to a 3d mesh. Because this mesh is 2d topologically but has 3d coordinates, this was not supported/broken. The functionality is now tested (first short "on_sphere" test) together with gmsh writing, parallel, checkpointing and flredecomp for this type of mesh.

Some more details of changes in Read_GMSH.F90:
- prevent reading of the gmsh file twice. This was because read_gmsh_simple() called identify_gmsh_file() to establish various dimensions, which it can only do by reading the entire file. Then an appropriate shape functon was allocated and read_gmsh_to_field was called which read the file again.
- only one version of read_gmsh_files to remain is the _simple version (as this was the only one in use). It
first reads in all the data, then establishes the necessary dimensions and then allocates the appropriate femtools
objects.
- the dimension of the read gmsh mesh is now established from the topology of the mesh, basically looking at the
highest dimension elements present, i.e. if there's tets or hexes it's 3d, otherwise 2d (1d gmsh reading is and remains untested). Previously this was done by looking at the z coordinate and seeing if it was small using some hard-coded definition of small=1e-9 (units?!).
- the dimension of the coordinate field that is returned (X%dim) is in most cases chosen to be the same as
the topological dimension. Only with the /geometry/spherical_earth option this is always 3 so that a 2d spherical
horizontal mesh (to be extruded inside fluidity) can be read. This was broken before. The other codes that use
embedded meshes (shallow water) use a specific option that only works for triangle, so that functionality should
remain unchanged.
- remove identify_gmsh_files()
- because of this remove identify_mesh_files() from the mesh_files module
- in all subroutines of the mesh_files module make format a required argument. Before, if the argument was not present a "guess" would be made. In some tools where no option tree is present (e.g. test_laplacian) this would always default to "triangle". These therefore now directly call read_triangle_files. In all other cases the format should be worked out. For output this can be a bit tricky. For checkpointing the same logic is followed as before (look at the option under the external mesh). For adaptivity debugging meshes it now always uses gmsh.
- identify_mesh_files() was only used in populate_state, in the case of inactive process (flredecomp from less
to more processes), where inactive processed needed to know the dimensions of the mesh. This is now communicated
via a MPI_Bcast()

[lp:~fluidity-core/fluidity/gmsh-on-sphere]

3955. By Stephan Kramer on 2012-03-28

Fix for silly unittest.

[Jon Hill (jon-hill)]

Code looks good at a quick eyeball.

Build queue is showing red:

unittests:
    Warnings: 1; tests = ['./test_empty_populate_state']
    Failures: 1; tests = ['./test_empty_populate_state']

Medium tests:
Summary of test problems with failures or warnings:
fladapt.xml: F

Any ideas?

[Lawrence Mitchell (wence)]

The fladapt issue is a real one, but it is generally applicable to all branches, rather than just this one I think. The error (with --enable-debugging) is something like:

Fortran runtime error: Actual string length is shorter than the declared one for dummy argument 'input' (-2123626932/7)

Here's a simplified test of the problematic code:

$ cat > foo.cpp <<\EOF
extern "C" {
#define F77_FUNC(name,NAME) name ## _
#define foo F77_FUNC(foo, FOO)
    void foo(const char *in, const int *len);
}

#include <iostream>
#include <stdlib.h>

using namespace std;
int main(int argc, char **argv)
{
    string in = "fladapt";
    int len = in.size();
    foo(in.c_str(), &len);
    return 0;
}
EOF

$ cat > foo.f90 <<\EOF
subroutine foo(input, length)
  integer, intent(in) :: length
  character(len = length), intent(in) :: input

  write(*,*)input
end subroutine foo
EOF

$ g++ -fbounds-check -O0 -g -c -o foo_main.o foo.cpp
$ gfortran -fbounds-check -O0 -g -c -o foo.o foo.f90
$ g++ -o foo -O0 -g foo.o foo_main.o -lgfortran
$ gdb ./foo
(gdb) r
Starting program: /home/lmitche4/tmp/foo
At line 1 of file foo.f90
Fortran runtime error: Actual string length is shorter than the declared one for dummy argument 'input' (-8200/7)

Program exited with code 02.
(gdb)

I think TRTD is either specified fixed length strings, or use iso_c_binding appropriately as is done in petsc_readnsolve. I think this is necessary for the majority of C++ main tools in tools/

[Jon Hill (jon-hill)]

I fixed a bunch of these a while back. Looks like there's a bunch more to fix...

Will add a bug report and set-up a branch

[Stephan Kramer (s-kramer)]

The unittest was already fixed in 3955. Just hasn't run again after that.

[Jon Hill (jon-hill)]

Just kicked the buildbot for this to check for greenness and will then approve.

[lp:~fluidity-core/fluidity/gmsh-on-sphere]

3956. By Jon Hill on 2012-05-01

Merge from trunk, only 1 conflict

[Stephan Kramer (s-kramer)]

Fladapt was fixed by Jon Hill in the trunk and merged in above. Will merge in again from trunk and restart the buildbot queue.

[lp:~fluidity-core/fluidity/gmsh-on-sphere]

3957. By Stephan Kramer on 2012-05-15

Merging in from trunk.

[Stephan Kramer (s-kramer)]

Ah, the fladapt failure actually didn't have anything to do with fortran/c binding issues, but was simply caused by my own stupidity. Fix is on the way.

[lp:~fluidity-core/fluidity/gmsh-on-sphere]

3958. By Stephan Kramer on 2012-05-16

Not a good idea to nullify a pointer just before trying to access its pointee.

[Jon Hill (jon-hill)]

Kicking buildbot again to test fix

[Jon Hill (jon-hill)]

Green. Approve.

review: Approve

=== modified file 'assemble/Adapt_State.F90'
--- assemble/Adapt_State.F90	2012-04-16 12:54:50 +0000
+++ assemble/Adapt_State.F90	2012-05-16 11:01:22 +0000
@@ -73,13 +73,12 @@
   use diagnostic_variables
   use diagnostic_fields_wrapper_new, only : calculate_diagnostic_variables_new => calculate_diagnostic_variables
   use pickers
+  use write_gmsh
 #ifdef HAVE_ZOLTAN
   use zoltan_integration
   use mpi_interfaces
 #endif
 
-  use mesh_files
-
   implicit none
 
   private
@@ -923,10 +922,10 @@
       output_positions => extract_vector_field(states(1), "Coordinate")
 
       if(isparallel()) then
-        call write_mesh_files(parallel_filename("first_timestep_adapted_mesh"), output_positions)
+        call write_gmsh_file(parallel_filename("first_timestep_adapted_mesh"), output_positions)
         call write_halos("first_timestep_adapted_mesh", output_positions%mesh)
       else
-        call write_mesh_files("first_timestep_adapted_mesh", output_positions)
+        call write_gmsh_file("first_timestep_adapted_mesh", output_positions)
       end if
 
     end if
@@ -1613,7 +1612,7 @@
       file_name = adapt_state_debug_file_name("adapted_mesh", mesh_dump_no)
       call find_mesh_to_adapt(states(1), mesh)
       positions = get_coordinate_field(states(1), mesh)
-      call write_mesh_files(file_name, positions)
+      call write_gmsh_file(file_name, positions)
       if(isparallel()) then
         file_name = adapt_state_debug_file_name("adapted_mesh", mesh_dump_no, add_parallel = .false.)  ! parallel extension is added by write_halos
         call write_halos(file_name, positions%mesh)
 
=== modified file 'assemble/Adapt_State_Prescribed.F90'
--- assemble/Adapt_State_Prescribed.F90	2011-07-22 15:34:14 +0000
+++ assemble/Adapt_State_Prescribed.F90	2012-05-16 11:01:22 +0000
@@ -109,8 +109,8 @@
     if(have_option(base_path // "/mesh/from_file")) then
       call get_option(base_path // "/mesh/from_file/format/name", format)
       call get_option("/geometry/quadrature/degree", quad_degree)
-          
-      new_positions = read_mesh_files(mesh_name, quad_degree = quad_degree)
+
+      new_positions = read_mesh_files(mesh_name, quad_degree = quad_degree, format = format)
     else
       ewrite(2, *) "Extracting new mesh from state"
     
 
=== modified file 'assemble/tests/test_empty_populate_state.F90'
--- assemble/tests/test_empty_populate_state.F90	2009-12-16 14:11:47 +0000
+++ assemble/tests/test_empty_populate_state.F90	2012-05-16 11:01:22 +0000
@@ -11,6 +11,7 @@
   logical :: fail
 
   is_active_process = .false.
+  no_active_processes = 0
   call load_options("data/empty-mesh.flml")
   call populate_state(states)
 
 
=== modified file 'femtools/Checkpoint.F90'
--- femtools/Checkpoint.F90	2011-07-19 15:52:43 +0000
+++ femtools/Checkpoint.F90	2012-05-16 11:01:22 +0000
@@ -1,4 +1,4 @@
-!    Copyright (C) 2006 Imperial College London and others.
+!    Copyrigh (C) 2006 Imperial College London and others.
 !
 !    Please see the AUTHORS file in the main source directory for a full list
 !    of copyright holders.
@@ -164,9 +164,9 @@
       & total_dete_lag
     type(vector_field), pointer :: vfield
 
-    integer(KIND=MPI_OFFSET_KIND) :: location_to_write, offset, offset_to_read
+    integer(KIND=MPI_OFFSET_KIND) :: location_to_write, offset
     integer, ALLOCATABLE, DIMENSION(:) :: status
-    integer :: nints, count, realsize, dimen, total_num_det, number_total_columns
+    integer :: nints, realsize, dimen, total_num_det, number_total_columns
     real, dimension(:), allocatable :: buffer
 
     ewrite(1, *) "Checkpointing detectors"
@@ -432,14 +432,12 @@
     logical, optional, intent(in) :: keep_initial_data
 
     type(vector_field), pointer:: position
-    character(len = FIELD_NAME_LEN) :: mesh_name
+    character(len = FIELD_NAME_LEN) :: mesh_name, mesh_format
     character(len = OPTION_PATH_LEN) :: mesh_path, mesh_filename
     integer :: i, n_meshes, stat1, stat2
-    type(mesh_type), pointer :: mesh
+    type(mesh_type), pointer :: mesh, external_mesh
     logical :: from_file, extruded
     
-    character(len = OPTION_PATH_LEN) :: currentMeshFormat
-
     assert(len_trim(prefix) > 0)
 
     n_meshes = option_count("/geometry/mesh")
@@ -467,12 +465,14 @@
         ! Update the options tree (required for options tree checkpointing)
         if (from_file) then
           call set_option_attribute(trim(mesh_path) // "/from_file/file_name", trim(mesh_filename))
+          call get_option(trim(mesh_path) // "/from_file/format/name", mesh_format)
         else if (extruded) then
 
-          ! Lunge forth with a wild stab at what the current mesh format is
-          call guess_external_mesh_format(currentMeshFormat)
+          ! the mesh format is determined from the external mesh
+          external_mesh => get_external_mesh(state)
+          call get_option(trim(external_mesh%option_path) // "/from_file/format/name", mesh_format)
 
-          call set_option_attribute(trim(mesh_path) // "/from_mesh/extrude/checkpoint_from_file/format/name", trim(currentMeshFormat), stat=stat1)
+          call set_option_attribute(trim(mesh_path) // "/from_mesh/extrude/checkpoint_from_file/format/name", trim(mesh_format), stat=stat1)
           call set_option_attribute(trim(mesh_path) // "/from_mesh/extrude/checkpoint_from_file/file_name", trim(mesh_filename), stat=stat2)
           if ((stat1/=SPUD_NO_ERROR .and. stat1/=SPUD_NEW_KEY_WARNING) .or. &
              & (stat2/=SPUD_NO_ERROR .and. stat2/=SPUD_NEW_KEY_WARNING)) then
@@ -487,13 +487,13 @@
           else
             position => extract_vector_field(state(1), trim(mesh%name)//"Coordinate")
           end if
-          call write_mesh_files(parallel_filename(mesh_filename), position)
+          call write_mesh_files(parallel_filename(mesh_filename), mesh_format, position)
           ! Write out the halos
           ewrite(2, *) "Checkpointing halos"
           call write_halos(mesh_filename, mesh)
         else
           ! Write out the mesh
-          call write_mesh_files(mesh_filename, state(1), mesh)
+          call write_mesh_files(mesh_filename, mesh_format, state(1), mesh)
         end if
      end if
    end do
 
=== modified file 'femtools/Mesh_Files.F90'
--- femtools/Mesh_Files.F90	2011-03-08 13:26:40 +0000
+++ femtools/Mesh_Files.F90	2012-05-16 11:01:22 +0000
@@ -62,8 +62,7 @@
   private
 
   interface read_mesh_files
-     module procedure read_mesh_files_to_field, &
-          read_mesh_files_to_state, read_mesh_simple
+     module procedure read_mesh_simple
   end interface
 
   interface write_mesh_files
@@ -71,8 +70,7 @@
           write_positions_to_file
   end interface
 
-  public :: read_mesh_files, identify_mesh_file, write_mesh_files
-  public :: guess_external_mesh_format
+  public :: read_mesh_files, write_mesh_files
 
 
 contains
@@ -82,143 +80,18 @@
   ! Read routines first
   ! --------------------------------------------------------------------------
 
-  subroutine identify_mesh_file(filename, dim, loc, nodes, elements, &
-       node_attributes, selements, selement_boundaries, format)
-    ! Discover the dimension and size of the mesh. Filename is 
-    ! the base name of the mesh file.
-    ! In parallel, filename must *include* the process number.
-
-    character(len=*), intent(in) :: filename
-    character(len=*), optional, intent(in) :: format
-    !! Dimension of mesh elements.
-    integer, intent(out), optional :: dim
-    !! Number of vertices of elements.
-    integer, intent(out), optional :: loc
-    !! Node and element counts.
-    integer, intent(out), optional :: nodes, elements
-    integer, intent(out), optional :: node_attributes
-    ! Surface element meta data
-    integer, optional, intent(out) :: selements
-    integer, optional, intent(out) :: selement_boundaries
-
-    character(len=OPTION_PATH_LEN) :: meshFormat
-
-    if( .not. present(format) ) then
-       call guess_external_mesh_format(meshFormat)
-    else
-       meshFormat = format
-    end if
-
-    ! Call appropriate subroutine depending upon format
-    select case( trim(meshFormat) )
-    case("triangle")
-       call identify_triangle_file(filename, dim, loc, nodes,&
-            elements, node_attributes, selements, selement_boundaries)
-
-    case("gmsh")
-       call identify_gmsh_file(filename, dim, loc, nodes,&
-            elements, node_attributes, selements, selement_boundaries)
-
-       ! Additional mesh format subroutines go here
-
-    case default
-       FLExit("Identifying mesh type "//format//" not supported within Fluidity")
-    end select
-
-  end subroutine identify_mesh_file
-
-
-  ! --------------------------------------------------------------------------
-
-
-  function read_mesh_files_to_field(filename, format, shape) result (field)
-    character(len=*), intent(in) :: filename
-    type(element_type), intent(in), target :: shape
-    type(vector_field)  :: field
-    character(len=OPTION_PATH_LEN) :: meshFormat
-
-    character(len=*), optional, intent(in) :: format
-
-    if( .not. present(format) ) then
-       call guess_external_mesh_format(meshFormat)
-    else
-       meshFormat = format
-    end if
-
-    select case( trim(meshFormat) )
-    case("triangle")
-       field = read_triangle_files(filename, shape)
-
-    case("gmsh")
-       field = read_gmsh_file(filename, shape)
-
-       ! Additional mesh format subroutines go here
-
-    case default
-       FLExit("Reading mesh type "//format//" not supported within Fluidity")
-    end select
-  end function read_mesh_files_to_field
-
-
-  ! --------------------------------------------------------------------------
-
-  function read_mesh_files_to_state(filename, shape, shape_type, &
-       n_states, format) &
-       result (result_state)
-
-    ! Filename is the base name of the mesh file without .ele, .msh, etc.
-    ! In parallel the filename must *not* include the process number.
-
-    character(len=*), intent(in) :: filename
-    type(element_type), intent(in), target :: shape
-    logical , intent(in):: shape_type
-    integer, intent(in), optional :: n_states
-
-    type(state_type)  :: result_state
-
-    character(len=*), optional, intent(in) :: format
-    character(len=option_path_len) :: meshFormat
-
-    if( .not. present(format) ) then
-       call guess_external_mesh_format(meshFormat)
-    else
-       meshFormat = format
-    end if
-
-
-    select case( trim(meshFormat) )
-    case("triangle")
-       result_state = read_triangle_files(filename, shape,shape_type,n_states)
-
-    case("gmsh")
-       result_state = read_gmsh_file(filename, shape,shape_type,n_states)
-
-       ! Additional mesh format subroutines go here
-
-    case default
-       FLExit("Reading mesh type "//format//" not supported within Fluidity")
-    end select
-
-  end function read_mesh_files_to_state
-
-
-  ! --------------------------------------------------------------------------
-
-  function read_mesh_simple(filename, quad_degree, &
-       quad_ngi, no_faces, &
-       quad_family, mdim, format ) &
+  function read_mesh_simple(filename, format, quad_degree, &
+       quad_ngi, quad_family, mdim) &
        result (field)
 
     ! A simpler mechanism for reading a mesh file into a field.
     ! In parallel the filename must *not* include the process number.
 
-    character(len=*), intent(in) :: filename
+    character(len=*), intent(in) :: filename, format
     ! The degree of the quadrature.
     integer, intent(in), optional, target :: quad_degree
     ! The degree of the quadrature.
     integer, intent(in), optional, target :: quad_ngi
-    ! Whether to add_faces on the resulting mesh.
-    logical, intent(in), optional :: no_faces
     ! What quadrature family to use
     integer, intent(in), optional :: quad_family
     ! Dimension of mesh
@@ -226,24 +99,17 @@
 
     type(vector_field) :: field
 
-    character(len=*), optional, intent(in) :: format
-    character(len=option_path_len) :: meshFormat
-
-    if( .not. present(format) ) then
-       call guess_external_mesh_format(meshFormat)
-    else
-       meshFormat = format
-    end if
-
-
-    select case( trim(meshFormat) )
+    select case( trim(format) )
     case("triangle")
-       field = read_triangle_files(filename, quad_degree, quad_ngi, &
-            no_faces, quad_family, mdim)
+       field = read_triangle_files(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, &
+            quad_family=quad_family, mdim=mdim)
 
     case("gmsh")
-       field = read_gmsh_file(filename, quad_degree, quad_ngi, &
-            no_faces, quad_family)
+       if (present(mdim)) then
+         FLExit("Cannot specify dimension for gmsh format")
+       end if
+       field = read_gmsh_file(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, &
+            quad_family=quad_family)
 
        ! Additional mesh format subroutines go here
 
@@ -260,22 +126,15 @@
   ! --------------------------------------------------------------------------
 
 
-  subroutine write_mesh_to_file(filename, state, mesh, format )
+  subroutine write_mesh_to_file(filename, format, state, mesh)
     ! Write out the supplied mesh to the specified filename as mesh files.
 
     character(len = *), intent(in) :: filename
-    character(len = *), intent(in), optional ::format
-    character(len = option_path_len) :: meshFormat
+    character(len = *), intent(in) :: format
     type(state_type), intent(in) :: state
     type(mesh_type), intent(in) :: mesh
 
-    if( .not. present(format) ) then
-       call guess_external_mesh_format(meshFormat)
-    else
-       meshFormat = format
-    end if
-
-    select case(meshFormat)
+    select case(format)
     case("triangle")
        call write_triangle_files(filename, state, mesh)
 
@@ -293,24 +152,13 @@
 
   ! --------------------------------------------------------------------------
 
-  subroutine write_positions_to_file(filename, positions, format)
+  subroutine write_positions_to_file(filename, format, positions)
     !!< Write out the mesh given by the position field in mesh files
     !!< In parallel, empty trailing processes are not written.
-    character(len=*), intent(in):: filename
+    character(len=*), intent(in):: filename, format
     type(vector_field), intent(in):: positions
-    !! "triangle" or "gmsh"
-    character(len=*), intent(in), optional :: format
-
-    character(len=option_path_len) :: meshFormat
-
-    if( .not. present(format) ) then
-       call guess_external_mesh_format(meshFormat)
-    else
-       meshFormat = format
-    end if
-
-
-    select case( trim(meshFormat) )
+
+    select case( trim(format) )
     case("triangle")
        call write_triangle_files( trim(filename), positions)
 
@@ -325,48 +173,5 @@
 
   end subroutine write_positions_to_file
 
-
-
-
-  ! --------------------------------------------------------------------------
-  ! Subroutine which finds external mesh format and puts it in 'meshFormat'
-  ! Follows similar logic to Field_Options::get_external_mesh()
-
-
-  subroutine guess_external_mesh_format(meshFormat)
-    character(len=*) :: meshFormat
-    character(len=OPTION_PATH_LEN) :: meshPath, formatPath
-    integer :: numMeshes, meshFound, i
-
-    numMeshes = option_count("/geometry/mesh")
-
-    ! Search for external meshes, and exit loop when found one
-    meshFound=0
-    do i = 0, numMeshes-1
-       meshPath = "/geometry/mesh["//int2str(i)//"]/from_file"
-
-       if(have_option(trim(meshPath))) then
-          meshFound=1
-          exit
-       end if
-    end do
-
-    ! We've found a mesh, now get its format
-    if (meshFound==1) then
-       formatPath = trim(meshPath) // "/format/name"
-
-       if(have_option(trim(formatPath))) then
-          call get_option( trim(formatPath), meshFormat )
-       else
-          ! If mesh exists but no format found, default to triangle
-          meshFormat = "triangle"
-       end if
-    else
-       ! If we can't find any external meshes...
-       meshFormat = "triangle"
-    end if
-
-  end subroutine guess_external_mesh_format
-
 end module mesh_files
 
 
=== modified file 'femtools/Read_GMSH.F90'
--- femtools/Read_GMSH.F90	2011-10-21 10:05:25 +0000
+++ femtools/Read_GMSH.F90	2012-05-16 11:01:22 +0000
@@ -44,160 +44,46 @@
   private
 
   interface read_gmsh_file
-     module procedure read_gmsh_file_to_field, read_gmsh_simple, &
-          read_gmsh_file_to_state
+     module procedure read_gmsh_simple
   end interface
 
-  public :: read_gmsh_file, identify_gmsh_file
+  public :: read_gmsh_file
+
+  integer, parameter:: GMSH_LINE=1, GMSH_TRIANGLE=2, GMSH_QUAD=3, GMSH_TET=4, GMSH_HEX=5, GMSH_NODE=15
 
 contains
 
   ! -----------------------------------------------------------------
-  ! GMSH version of triangle equivalent.
-
-  subroutine identify_gmsh_file(filename, numDimenOut, locOut, &
-       numNodesOut, numElementsOut, &
-       nodeAttributesOut, selementsOut, boundaryFlagOut)
-    ! Discover the dimension and size of the GMSH mesh.
-    ! Filename is the base name of the file without .msh.
-    ! In parallel, filename must *include* the process number.
-
-    character(len=*), intent(in) :: filename
-    character(len=option_path_len) :: lfilename
-
-    !! Number of vertices of elements.
-    integer, intent(out), optional :: numDimenOut, locOut, numNodesOut
-    integer, intent(out), optional :: numElementsOut, nodeAttributesOut
-    integer, intent(out), optional :: selementsOut, boundaryFlagOut
-
-
-    logical :: fileExists
-
-    integer :: fd, gmshFormat
-    type(GMSHnode), pointer :: nodes(:)
-    type(GMSHelement), pointer :: elements(:), faces(:)
-    integer :: numElements, boundaryFlag, numNodes, numDimen
-    integer :: loc, effDimen, nodeAttributes
-    integer :: i
-
-    logical :: haveBounds, haveInternalBounds
-
-
-    lfilename = trim(filename) // ".msh"
-
-    ! Read node file header
-    inquire(file = trim(lfilename), exist = fileExists)
-    if(.not. fileExists) then
-       FLExit("gmsh file " // trim(lfilename) // " not found")
-    end if
-
-    ewrite(2, *) "Opening " // trim(lfilename) // " for reading."
-    fd = free_unit()
-    open(unit = fd, file = trim(lfilename), &
-         err=42, action="read", access="stream", form="formatted" )
-
-    call read_header(fd, lfilename, gmshFormat)
-
-    ! Read in the nodes
-    call read_nodes_coords( fd, lfilename, gmshFormat, nodes )
-
-    numDimen = mesh_dimensions( nodes )
-    numNodes = size(nodes)
-
-    ! Read in elements
-    call read_faces_and_elements( fd, lfilename, gmshFormat, &
-         elements, faces )
-
-
-    ! Try reading in node column ID data (if there is any)
-    call read_node_column_IDs( fd, lfilename, gmshFormat, nodes)
-
-    close( fd )
-
-    ! We're assuming all elements have the same number of vertices/nodes
-    loc = size(elements(1)%nodeIDs)
-
-    do i=1, numNodes
-       if( nodes(i)%columnID>0 ) nodeAttributes=1
-    end do
-
-    ! NOTE:  'boundaries' variable
-    ! (see Read_Triangle.F90) is a flag which indicates whether
-    ! faces have boundaries (physical IDs). Values:
-    ! =0  : no boundaries
-    ! =1  : boundaries, normal
-    ! =2 : boundaries, periodic mesh (internal boundary)
-
-    haveBounds=.false.
-    haveInternalBounds=.false.
-
-    do i=1, size(faces)
-       if(faces(i)%numTags > 0) then
-          ! We have boundaries,
-          haveBounds=.true.
-
-          ! We have internal boundaries (at the join of a period mesh)
-          if(faces(i)%numTags >= 4) then
-             if (faces(i)%tags(4) > 0) then
-                haveInternalBounds=.true.
-             end if
-          end if
-       end if
-    end do
-
-    if(haveInternalBounds) then
-       boundaryFlag=2
-    elseif(haveBounds) then
-       boundaryFlag=1
-    else
-       boundaryFlag=0
-    end if
-
-    if( numDimen.eq.2 .and. have_option("/geometry/spherical_earth/") ) then
-       effDimen = numDimen+1
-    else
-       effDimen = numDimen
-    end if
-
-
-    ! Return optional variables requested
-
-    if(present(nodeAttributesOut)) nodeAttributesOut=nodeAttributes
-    if(present(numDimenOut)) numDimenOut=effDimen
-    if(present(numElementsOut)) numElementsOut=numElements
-    if(present(locOut)) locOut=loc
-    if(present(boundaryFlagOut)) boundaryFlagOut=boundaryFlag
-
-    deallocate(nodes)
-
-    return
-
-42  FLExit("Unable to open "//trim(lfilename))
-
-  end subroutine identify_gmsh_file
-
-
-
-  ! -----------------------------------------------------------------
   ! The main function for reading GMSH files
 
-  function read_gmsh_file_to_field(filename, shape) result (field)
-    ! Filename is the base name of the GMSH file without .msh 
-    ! In parallel the filename must *not* include the process number.
+  function read_gmsh_simple( filename, quad_degree, &
+       quad_ngi, quad_family ) &
+       result (field)
+    !!< Read a GMSH file into a coordinate field.
+    !!< In parallel the filename must *not* include the process number.
 
     character(len=*), intent(in) :: filename
-    type(element_type), intent(in), target :: shape
-    type(vector_field)  :: field
+    !! The degree of the quadrature.
+    integer, intent(in), optional, target :: quad_degree
+    !! The degree of the quadrature.
+    integer, intent(in), optional, target :: quad_ngi
+    !! What quadrature family to use
+    integer, intent(in), optional :: quad_family
+    !! result: a coordinate field
+    type(vector_field) :: field
+
+    type(quadrature_type):: quad
+    type(element_type):: shape
+    type(mesh_type):: mesh
 
     integer :: fd
     integer,  pointer, dimension(:) :: sndglno, boundaryIDs, faceOwner
 
     character(len = parallel_filename_len(filename)) :: lfilename
     integer :: loc, sloc
-    type(mesh_type) :: mesh
     integer :: numNodes, numElements, numFaces
     logical :: haveBounds, haveElementOwners, haveRegionIDs
-    integer :: numDimen, effDimen
+    integer :: dim, coordinate_dim
     integer :: gmshFormat
     integer :: n, d, e, f, nodeID
 
@@ -225,11 +111,9 @@
     ! Read in the nodes
     call read_nodes_coords( fd, lfilename, gmshFormat, nodes )
 
-    numDimen = mesh_dimensions( nodes )
-
     ! Read in elements
     call read_faces_and_elements( fd, lfilename, gmshFormat, &
-         elements, faces )
+         elements, faces, dim)
 
     call read_node_column_IDs( fd, lfilename, gmshFormat, nodes )
 
@@ -285,39 +169,51 @@
     
     end if
     
-    if( numDimen.eq.2 .and. have_option("/geometry/spherical_earth/") ) then
-       effDimen = numDimen+1
-    else
-       effDimen = numDimen
-    end if
-
-
+    if(have_option("/geometry/spherical_earth/") ) then
+      ! on the sphere the input mesh may be 2d (extrusion), or 3d but
+      ! Coordinate is always 3-dimensional
+      coordinate_dim  = 3
+    else
+      coordinate_dim  = dim
+    end if
+
+    loc = size( elements(1)%nodeIDs )
+    if (numFaces>0) then
+      sloc = size( faces(1)%nodeIDs )
+    else
+      sloc = 0
+    end if
+
+    ! Now construct within Fluidity data structures
+
+    if (present(quad_degree)) then
+       quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family)
+    else if (present(quad_ngi)) then
+       quad = make_quadrature(loc, dim, ngi=quad_ngi, family=quad_family)
+    else
+       FLAbort("Need to specify either quadrature degree or ngi")
+    end if
+    shape=make_element_shape(loc, dim, 1, quad)
     call allocate(mesh, numNodes, numElements, shape, name="CoordinateMesh")
-    call allocate( field, effDimen, mesh, name="Coordinate")
-    call deallocate( mesh )
-
-    ! Now construct within Fluidity data structures
+    call allocate( field, coordinate_dim, mesh, name="Coordinate")
+
+    ! deallocate our references of mesh, shape and quadrature:
+    call deallocate(mesh)
+    call deallocate(shape)
+    call deallocate(quad)
+
 
     if (haveRegionIDs) then
       allocate( field%mesh%region_ids(numElements) )
     end if
-    if(nodes(1)%columnID.ge.0)  allocate(field%mesh%columns(1:numNodes))
-
-    loc = size( elements(1)%nodeIDs )
-    if (numFaces>0) then
-      sloc = size( faces(1)%nodeIDs )
-    else
-      sloc = 0
-    end if
-
-    assert(loc==shape%loc)
+    if(nodes(1)%columnID>=0)  allocate(field%mesh%columns(1:numNodes))
 
     ! Loop round nodes copying across coords and column IDs to field mesh,
     ! if they exist
     do n=1, numNodes
 
        nodeID = nodes(n)%nodeID
-       forall (d = 1:effDimen)
+       forall (d = 1:field%dim)
           field%val(d,nodeID) = nodes(n)%x(d)
        end forall
 
@@ -380,86 +276,8 @@
 
 43  FLExit("Unable to open "//trim(lfilename))
 
-  end function read_gmsh_file_to_field
-
-
-
-
-
-
-  ! -----------------------------------------------------------------
-  ! Simplified interface to reading GMSH files
-  function read_gmsh_simple( filename, quad_degree, &
-       quad_ngi, no_faces, quad_family ) &
-       result (field)
-    !!< A simpler mechanism for reading a GMSH file into a field.
-    !!< In parallel the filename must *not* include the process number.
-
-    character(len=*), intent(in) :: filename
-    !! The degree of the quadrature.
-    integer, intent(in), optional, target :: quad_degree
-    !! The degree of the quadrature.
-    integer, intent(in), optional, target :: quad_ngi
-    !! Whether to add_faces on the resulting mesh.
-    logical, intent(in), optional :: no_faces
-    !! What quadrature family to use
-    integer, intent(in), optional :: quad_family
-
-    type(vector_field) :: field
-    type(quadrature_type) :: quad
-    type(element_type) :: shape
-
-    integer :: dim, loc
-
-
-    if(isparallel()) then
-       call identify_gmsh_file(parallel_filename(filename), dim, loc)
-    else
-       call identify_gmsh_file(filename, dim, loc)
-    end if
-
-
-    if (present(quad_degree)) then
-       quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family)
-    else if (present(quad_ngi)) then
-       quad = make_quadrature(loc, dim, ngi=quad_ngi, family=quad_family)
-    else
-       FLAbort("Need to specify either quadrature degree or ngi")
-    end if
-
-
-    shape=make_element_shape(loc, dim, 1, quad)
-
-    field=read_gmsh_file(filename, shape)
-
-    ! deallocate our references of shape and quadrature:
-    call deallocate_element(shape)
-    call deallocate(quad)
-
   end function read_gmsh_simple
 
-
-
-
-  ! -----------------------------------------------------------------
-  ! Read GMSH file to state object.
-
-  function read_gmsh_file_to_state(filename, shape,shape_type,n_states) &
-       result (result_state)
-    ! Filename is the base name of the GMSH file without .node or .ele.
-    ! In parallel the filename must *not* include the process number.
-
-    character(len=*), intent(in) :: filename
-    type(element_type), intent(in), target :: shape
-    logical , intent(in):: shape_type
-    integer, intent(in), optional :: n_states
-    type(state_type)  :: result_state
-
-    FLAbort("read_gmsh_file_to_state() not implemented yet")
-
-  end function read_gmsh_file_to_state
-
-
   ! -----------------------------------------------------------------
   ! Read through the head to decide whether binary or ASCII, and decide
   ! whether this looks like a GMSH mesh file or not.
@@ -666,37 +484,17 @@
 
 
   ! -----------------------------------------------------------------
-  ! Guesstimate mesh dimensions. If any node z-coords are non-zero, then
-  ! it's a 3D mesh; otherwise 2D. Not pretty, but there seems no other
-  ! way to do this.
-
-  integer function mesh_dimensions( nodes )
-    type(GMSHnode), pointer :: nodes(:)
-    real, pointer :: absZ(:)
-
-    mesh_dimensions = 3
-
-    allocate ( absZ(size(nodes)) )
-    absZ(:) = nodes(:)%x(3)
-
-    if( all( absZ > -verySmall .and. absZ < verySmall) ) mesh_dimensions=2
-
-    deallocate(absZ)
-  end function mesh_dimensions
-
-
-
-
-  ! -----------------------------------------------------------------
-  ! Read in element header data
+  ! Read in element header data and establish topological dimension
 
   subroutine read_faces_and_elements( fd, filename, gmshFormat, &
-       elements, faces )
-
-    integer :: fd, gmshFormat
-    character(len=*) :: filename
-
-    type(GMSHelement), pointer :: allElements(:), elements(:), faces(:)
+       elements, faces, dim)
+
+    integer, intent(in) :: fd, gmshFormat
+    character(len=*), intent(in) :: filename
+    type(GMSHelement), pointer :: elements(:), faces(:)
+    integer, intent(out) :: dim
+
+    type(GMSHelement), pointer :: allElements(:)
 
     integer :: numAllElements
     character(len=longStringLen) :: charBuf
@@ -708,14 +506,14 @@
 
 
     read(fd,*) charBuf
-    if( trim(charBuf) .ne. "$Elements" ) then
+    if( trim(charBuf)/="$Elements" ) then
        FLExit("Error: cannot find '$Elements' in GMSH mesh file")
     end if
 
     read(fd,*) numAllElements
 
     ! Sanity check.
-    if(numAllElements .lt. 1) then
+    if(numAllElements<1) then
        FLExit("Error: number of elements in GMSH file < 1")
     end if
 
@@ -799,23 +597,21 @@
             allElements(e)%type )
 
        select case ( allElements(e)%type )
-          ! A line
-       case (1)
+       case (GMSH_LINE)
           numEdges = numEdges+1
-          ! A triangle
-       case (2)
+       case (GMSH_TRIANGLE)
           numTriangles = numTriangles+1
-          ! A quad
-       case (3)
+       case (GMSH_QUAD)
           numQuads = numQuads+1
-       case (4)
+       case (GMSH_TET)
           numTets = numTets+1
-       case (5)
+       case (GMSH_HEX)
           numHexes = numHexes+1
-       case (15)
+       case (GMSH_NODE)
           ! Do nothing
        case default
-          FLExit("read_faces_and_elements(): unsupported element type")
+          ewrite(0,*) "element id,type: ", allElements(e)%elementID, allElements(e)%type
+          FLExit("Unsupported element type in gmsh .msh file")
        end select
 
     end do
@@ -833,29 +629,40 @@
     ! meshes are verboten:
     !   tet/hex, tet/quad, triangle/hex and triangle/quad
 
-    if (numTets .gt. 0) then
+    if (numTets>0) then
        numElements = numTets
-       elementType = 4
+       elementType = GMSH_TET
        numFaces = numTriangles
-       faceType = 2
+       faceType = GMSH_TRIANGLE
+       dim = 3
+       if (numQuads>0 .or. numHexes>0) then
+         FLExit("Cannot combine hexes or quads with tetrahedrals in one gmsh .msh file")
+       end if
 
-    elseif (numTriangles .gt. 0) then
+    elseif (numTriangles>0) then
        numElements = numTriangles
-       elementType = 2
+       elementType = GMSH_TRIANGLE
        numFaces = numEdges
-       faceType = 1
+       faceType = GMSH_LINE
+       dim = 2
+       if (numQuads>0 .or. numHexes>0) then
+         FLExit("Cannot combine hexes or quads with triangles in one gmsh .msh file")
+       end if
 
     elseif (numHexes .gt. 0) then
        numElements = numHexes
-       elementType = 5
+       elementType = GMSH_HEX
        numFaces = numQuads
-       faceType = 3
+       faceType = GMSH_QUAD
+       dim = 3
 
     elseif (numQuads .gt. 0) then
        numElements = numQuads
-       elementType = 3
+       elementType = GMSH_QUAD
        numFaces = numEdges
-       faceType = 1
+       faceType = GMSH_LINE
+       dim = 2
+
     else
        FLExit("Unsupported mixture of face/element types")
     end if
 
=== modified file 'femtools/Read_Triangle.F90'
--- femtools/Read_Triangle.F90	2011-03-08 13:26:40 +0000
+++ femtools/Read_Triangle.F90	2012-05-16 11:01:22 +0000
@@ -42,8 +42,7 @@
   private
 
   interface read_triangle_files
-     module procedure read_triangle_files_to_field, &
-          read_triangle_files_to_state, read_triangle_simple
+     module procedure read_triangle_files_to_field, read_triangle_simple
   end interface
 
   public :: read_triangle_files, identify_triangle_file, read_elemental_mappings, read_triangle_serial
@@ -414,167 +413,6 @@
     
   end function read_triangle_files_to_field
 
-  function read_triangle_files_to_state(filename, shape,shape_type,n_states) result (result_state)
-    !!< Filename is the base name of the triangle file without .node or .ele.
-    !!< In parallel the filename must *not* include the process number.
-    
-    character(len=*), intent(in) :: filename
-    type(element_type), intent(in), target :: shape
-    logical , intent(in):: shape_type
-    integer, intent(in), optional :: n_states
-    type(state_type)  :: result_state
-    type(vector_field) :: field
-    type(scalar_field), allocatable, dimension(:) :: attribs
-    type(scalar_field) :: bndry_mark
-
-    character(len = parallel_filename_len(filename)) :: lfilename
-    integer :: node_unit, ele_unit
-    real, allocatable, dimension(:) :: read_buffer
-    character(len=6) :: write_buffer
-    integer :: i, j, nodes, dim, xdim, node_attributes, boundaries,&
-         & ele_attributes, loc, elements
-    integer, allocatable, dimension(:) :: node_order
-
-    type(mesh_type), dimension(:), allocatable, save :: meshes
-    integer, save :: counter=0
-
-    assert(shape_type)
-    
-    ! If running in parallel, add the process number
-    if(isparallel()) then
-      lfilename = parallel_filename(filename)
-    else
-      lfilename = trim(filename)
-    end if
-
-    node_unit=free_unit()
-
-    counter=counter+1
-    if (counter==1) then
-       if (present(n_states)) then
-          allocate(meshes(n_states))
-       else
-          allocate(meshes(1))
-       end if
-    end if
-    assert(counter .le. size(meshes))
-
-    call nullify(result_state)
-    call deallocate(result_state)
-
-    ewrite(2, *) "Opening "//trim(lfilename)//".node for reading."
-    ! Open node file
-    open(unit=node_unit, file=trim(lfilename)//".node", err=42, action="read")
-    
-    ! Read node file header.
-    read (node_unit, *) nodes, xdim, node_attributes, boundaries
-
-    ele_unit=free_unit()
-
-    ewrite(2, *) "Opening "//trim(lfilename)//".ele for reading."
-    ! Open element file
-    open(unit=ele_unit, file=trim(lfilename)//".ele", err=43, action="read")
-    
-    ! Read element file header.
-    read (ele_unit, *) elements, loc, ele_attributes
-
-    assert(loc==shape%loc)
-    allocate(node_order(loc))
-    select case(loc)
-    case(3)
-       node_order = (/1,2,3/)
-    case(6)
-       node_order = (/1,6,2,5,4,3/)
-    case default
-       do j=1,loc
-          node_order(:)=j
-       enddo
-    end select
-    
-    call allocate(meshes(counter), nodes, elements, shape, name="CoordinateMesh")
-    call allocate(field, xdim, meshes(counter), name=filename//' Coordinate')
-    allocate(attribs(node_attributes))
-    do j=1,node_attributes
-       write(write_buffer,'(i0)') j
-       call allocate(attribs(j),meshes(counter),&
-            name='Attribute '//trim(write_buffer))
-    end do
-       call allocate(bndry_mark,meshes(counter),&
-            name='Boundary Marker')
-
-    allocate(read_buffer(xdim+node_attributes+boundaries+1))
-
-    do i=1,nodes
-       read(node_unit,*) read_buffer
-
-       forall (j=1:xdim)
-          field%val(j,i)=read_buffer(j+1)
-       end forall
-
-       forall (j=1:node_attributes)
-          attribs(j)%val(i)=read_buffer(j+1+xdim)
-       end forall
-
-       if (boundaries>0) then
-          bndry_mark%val(i)=read_buffer(1+1+xdim+node_attributes)
-       end if
-
-    end do
-
-
-    deallocate(read_buffer)
-    allocate(read_buffer(loc+ele_attributes+1))
-
-    if(ele_attributes==1) then  ! this assumes the element attribute is a region id
-      allocate(field%mesh%region_ids(1:elements))
-      do j=1,node_attributes
-        allocate(attribs(j)%mesh%region_ids(1:elements))
-      end do
-      allocate(bndry_mark%mesh%region_ids(1:elements))
-      allocate(meshes(counter)%region_ids(1:elements))
-    end if
-
-    do i=1,elements
-       read(ele_unit,*) read_buffer
-
-
-       field%mesh%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1)
-       forall (j=1:node_attributes)
-          attribs(j)%mesh%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1)
-       end forall
-       bndry_mark%mesh%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1)
-       meshes(counter)%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1)
-
-       if(ele_attributes==1) then
-          field%mesh%region_ids(i)=read_buffer(loc+2)
-          forall (j=1:node_attributes)
-              attribs(j)%mesh%region_ids(i)=read_buffer(loc+2)
-          end forall
-          bndry_mark%mesh%region_ids(i)=read_buffer(loc+2)
-          meshes(counter)%region_ids(i)=read_buffer(loc+2)
-       end if
-    end do
-    
-    close(node_unit)
-    close(ele_unit)
-
-
-    call insert(result_state,field,"Coordinate")    
-    do j=1,node_attributes
-       write(write_buffer,'(i0)') j
-       call insert(result_state,attribs(j),"Attributes "//trim(write_buffer))
-    end do
-    if(boundaries>0) call insert(result_state,bndry_mark,"Boundary Marker")
-    call insert(result_state,meshes(counter),filename)
-
-    return
-
-42  FLExit("Unable to open "//trim(lfilename)//".node")
-
-43  FLExit("Unable to open "//trim(lfilename)//".ele")
-
-  end function read_triangle_files_to_state
-
   function read_triangle_simple(filename, quad_degree, quad_ngi, no_faces, quad_family, mdim) result (field)
     !!< A simpler mechanism for reading a triangle file into a field.
     !!< In parallel the filename must *not* include the process number.
 
=== modified file 'preprocessor/Populate_State.F90'
--- preprocessor/Populate_State.F90	2012-01-30 15:36:17 +0000
+++ preprocessor/Populate_State.F90	2012-05-16 11:01:22 +0000
@@ -174,6 +174,7 @@
     character(len=OPTION_PATH_LEN) :: mesh_path, mesh_file_name,&
          & mesh_file_format, from_file_path
     integer, dimension(:), pointer :: coplanar_ids
+    integer, dimension(3) :: mesh_dims
     integer :: i, j, nmeshes, nstates, quad_degree, stat
     type(element_type), pointer :: shape
     type(quadrature_type), pointer :: quad
@@ -219,22 +220,96 @@
           call get_option("/geometry/quadrature/degree", quad_degree)
           quad_family = get_quad_family()
 
+
+          if (is_active_process) then
+            select case (mesh_file_format)
+            case ("triangle", "gmsh")
+              ! Get mesh dimension if present
+              call get_option(trim(mesh_path)//"/from_file/dimension", mdim, stat)
+              ! Read mesh
+              if(stat==0) then
+                 position=read_mesh_files(trim(mesh_file_name), &
+                      quad_degree=quad_degree, &
+                      quad_family=quad_family, mdim=mdim, &
+                      format=mesh_file_format)
+              else
+                 position=read_mesh_files(trim(mesh_file_name), &
+                      quad_degree=quad_degree, &
+                      quad_family=quad_family, &
+                      format=mesh_file_format)
+              end if
+              mesh=position%mesh
+            case ("vtu")
+              position_ptr => vtk_cache_read_positions_field(mesh_file_name)
+              ! No hybrid mesh support here
+              assert(ele_count(position_ptr) > 0)
+              dim = position_ptr%dim
+              loc = ele_loc(position_ptr, 1)
+
+              ! Generate a copy, and swap the quadrature degree
+              ! Note: Even if positions_ptr has the correct quadrature degree, it
+              ! won't have any faces and hence a copy is still required (as
+              ! add_faces is a construction routine only)
+              allocate(quad)
+              allocate(shape)
+              quad = make_quadrature(loc, dim, degree = quad_degree, family=quad_family)
+              shape = make_element_shape(loc, dim, 1, quad)
+              call allocate(mesh, nodes = node_count(position_ptr), elements = ele_count(position_ptr), shape = shape, name = position_ptr%mesh%name)
+              do j = 1, ele_count(mesh)
+                 call set_ele_nodes(mesh, j, ele_nodes(position_ptr%mesh, j))
+              end do
+              call add_faces(mesh)
+              call allocate(position, dim, mesh, position_ptr%name)
+              call set(position, position_ptr)
+              call deallocate(mesh)
+              call deallocate(shape)
+              call deallocate(quad)
+              deallocate(quad)
+              deallocate(shape)
+
+              mesh = position%mesh
+            case default
+              ewrite(-1,*) trim(mesh_file_format), " is not a valid format for a mesh file"
+              FLAbort("Invalid format for mesh file")
+            end select
+         end if
+
+          if (no_active_processes /= getnprocs()) then
+            ! not all processes are active, they need to be told the mesh dimensions
+
+            ! receive the mesh dimension from rank 0
+            if (getrank()==0) then
+              if (is_active_process) then
+                ! normally rank 0 should always be active, so it knows the dimensions
+                mesh_dims(1)=mesh_dim(mesh)
+                mesh_dims(2)=ele_loc(mesh,1)
+                if (associated(mesh%columns)) then
+                  mesh_dims(3)=1
+                else
+                  mesh_dims(3)=0
+                end if
+              else
+                ! this is a special case for a unit test with 1 inactive process
+                call get_option('/geometry/dimension', mesh_dims(1))
+                mesh_dims(2)=mesh_dims(1)+1
+                mesh_dims(3)=0
+              end if
+            end if
+            call MPI_bcast(mesh_dims, 3, getpinteger(), 0, MPI_COMM_FEMTOOLS, stat)
+          end if
+
+          
           if (.not. is_active_process) then
             ! is_active_process records whether we have data on disk or not
             ! see the comment in Global_Parameters. In this block, 
             ! we want to allocate an empty mesh and positions.
 
+            dim=mesh_dims(1)
+            loc=mesh_dims(2)
+            column_ids=mesh_dims(3)
+
             allocate(quad)
             allocate(shape)
-            if (no_active_processes == 1) then
-              call identify_mesh_file(trim(mesh_file_name), dim, loc, &
-                node_attributes=column_ids, &
-                format=mesh_file_format )
-            else
-              call identify_mesh_file(trim(mesh_file_name) // "_0", dim, loc, &
-                node_attributes=column_ids, &
-                format=mesh_file_format )
-            end if
             quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family)
             shape=make_element_shape(loc, dim, 1, quad)
             call allocate(mesh, nodes=0, elements=0, shape=shape, name="EmptyMesh")
@@ -252,58 +327,8 @@
             
             deallocate(quad)
             deallocate(shape)
-
-         else if(trim(mesh_file_format)=="triangle" &
-              .or. trim(mesh_file_format)=="gmsh") then
-            ! Get mesh dimension if present
-            call get_option(trim(mesh_path)//"/from_file/dimension", mdim, stat)
-            ! Read mesh
-            if(stat==0) then
-               position=read_mesh_files(trim(mesh_file_name), &
-                    quad_degree=quad_degree, &
-                    quad_family=quad_family, mdim=mdim, &
-                    format=mesh_file_format)
-            else
-               position=read_mesh_files(trim(mesh_file_name), &
-                    quad_degree=quad_degree, &
-                    quad_family=quad_family, &
-                    format=mesh_file_format)
-            end if
-            mesh=position%mesh
-         else if(trim(mesh_file_format) == "vtu") then
-            position_ptr => vtk_cache_read_positions_field(mesh_file_name)
-            ! No hybrid mesh support here
-            assert(ele_count(position_ptr) > 0)
-            dim = position_ptr%dim
-            loc = ele_loc(position_ptr, 1)
-
-            ! Generate a copy, and swap the quadrature degree
-            ! Note: Even if positions_ptr has the correct quadrature degree, it
-            ! won't have any faces and hence a copy is still required (as
-            ! add_faces is a construction routine only)
-            allocate(quad)
-            allocate(shape)
-            quad = make_quadrature(loc, dim, degree = quad_degree, family=quad_family)
-            shape = make_element_shape(loc, dim, 1, quad)
-            call allocate(mesh, nodes = node_count(position_ptr), elements = ele_count(position_ptr), shape = shape, name = position_ptr%mesh%name)
-            do j = 1, ele_count(mesh)
-               call set_ele_nodes(mesh, j, ele_nodes(position_ptr%mesh, j))
-            end do
-            call add_faces(mesh)
-            call allocate(position, dim, mesh, position_ptr%name)
-            call set(position, position_ptr)
-            call deallocate(mesh)
-            call deallocate(shape)
-            call deallocate(quad)
-            deallocate(quad)
-            deallocate(shape)
-
-            mesh = position%mesh
-         else
-            ewrite(-1,*) trim(mesh_file_format), " is not a valid format for a mesh file"
-            FLAbort("Invalid format for mesh file")
-         end if
-          
+          end if
+
           ! if there is a derived mesh which specifies periodic bcs 
           ! to be *removed*, we assume the external mesh is periodic
           mesh%periodic = option_count("/geometry/mesh/from_mesh/&
 
=== added directory 'tests/spherical_patch'
=== added file 'tests/spherical_patch/Makefile'
--- tests/spherical_patch/Makefile	1970-01-01 00:00:00 +0000
+++ tests/spherical_patch/Makefile	2012-05-16 11:01:22 +0000
@@ -0,0 +1,11 @@
+GMSH_GEO_FILE=src/spherical_segment.geo
+NPROCS=3
+input: clean
+	gmsh -2 -bin $(GMSH_GEO_FILE) -o patch.msh
+	../../bin/fldecomp -n $(NPROCS) -s -m gmsh patch
+
+clean:
+	rm -rf *.msh *.halo *.vtu *.stat fluidity.*
+	rm -rf spherical_segment_* *.pvtu
+	rm -rf flredecomp.* matrixdump*
+
 
=== added file 'tests/spherical_patch/spherical_patch.xml'
--- tests/spherical_patch/spherical_patch.xml	1970-01-01 00:00:00 +0000
+++ tests/spherical_patch/spherical_patch.xml	2012-05-16 11:01:22 +0000
@@ -0,0 +1,47 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>spherical patch<comment>Tests gmsh reading and writing on the sphere with 2d horizontal spherical input mesh, checkpointing and flredecomp.</comment></name>
+  <owner userid="skramer"/>
+  <problem_definition length="short" nprocs="1">
+    <command_line>mpirun -np 3 fluidity -v2 -l spherical_segment.flml &amp;&amp;
+mpirun -np 3 ../../bin/flredecomp -v -l -i 3 -o 1 spherical_segment_2_checkpoint serial &amp;&amp;
+fluidity -v2 -l serial.flml</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="volume_before" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('spherical_segment.stat')
+# this only works because we don't move the free surface
+# we actually compute H \int_top \eta, i.e. the volume difference with
+# equilibrium state multiplied by depth H
+volume_before=stat['Fields']['FreeSurface']['integral']</variable>
+    <variable name="volume_after" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('spherical_segment_checkpoint.stat')
+# this only works because we don't move the free surface
+# we actually compute H \int_top \eta, i.e. the volume difference with
+# equilibrium state multiplied by depth H
+volume_after=stat['Fields']['FreeSurface']['integral']</variable>
+    <variable name="solvers_converged" language="python">import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files</variable>
+    <variable name="max_velocity_before" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('spherical_segment.stat')
+max_velocity_before=stat['Fields']['Velocity%magnitude']['max']</variable>
+    <variable name="max_velocity_after" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('spherical_segment_checkpoint.stat')
+max_velocity_after=stat['Fields']['Velocity%magnitude']['max']</variable>
+  </variables>
+  <pass_tests>
+    <test name="volume_conservation_run1" language="python">from fluidity_tools import compare_variable
+compare_variable(volume_before[0], volume_before[-1], 1e-5)</test>
+    <test name="volume_conservation_run12" language="python">from fluidity_tools import compare_variable
+# compares volume at the end of first run with that of end first timestep second run
+compare_variable(volume_before[-1], volume_after[0], 1e-5)</test>
+    <test name="volume_conservation_run2" language="python">from fluidity_tools import compare_variable
+compare_variable(volume_after[0], volume_after[-1], 1e-5)</test>
+    <test name="solvers_converged" language="python">assert solvers_converged</test>
+    <test name="max_velocity_run1" language="python"># make sure we're actually doing something and not blowing up
+assert all(max_velocity_before&gt;0.01) and all(max_velocity_before&lt;0.02)</test>
+    <test name="max_velocity_run2" language="python"># make sure we're actually doing something and not blowing up
+assert all(max_velocity_after&gt;0.01) and all(max_velocity_after&lt;0.02)</test>
+  </pass_tests>
+</testproblem>
 
=== added file 'tests/spherical_patch/spherical_segment.flml'
--- tests/spherical_patch/spherical_segment.flml	1970-01-01 00:00:00 +0000
+++ tests/spherical_patch/spherical_segment.flml	2012-05-16 11:01:22 +0000
@@ -0,0 +1,502 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">spherical_segment</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">oceans</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_mesh>
+        <mesh name="BaseMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="BaseMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="BaseMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="InputMesh">
+      <from_file file_name="patch">
+        <format name="gmsh"/>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="BaseMesh">
+      <from_mesh>
+        <mesh name="InputMesh"/>
+        <extrude>
+          <regions name="WholeMesh">
+            <bottom_depth>
+              <constant>
+                <real_value rank="0">90</real_value>
+              </constant>
+            </bottom_depth>
+            <sizing_function>
+              <constant>
+                <real_value rank="0">30</real_value>
+              </constant>
+            </sizing_function>
+            <top_surface_id>
+              <integer_value rank="0">1</integer_value>
+            </top_surface_id>
+            <bottom_surface_id>
+              <integer_value rank="0">2</integer_value>
+            </bottom_surface_id>
+          </regions>
+        </extrude>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="P0Mesh">
+      <from_mesh>
+        <mesh name="BaseMesh"/>
+        <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>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+      <surface_degree>
+        <integer_value rank="0">3</integer_value>
+      </surface_degree>
+    </quadrature>
+    <spherical_earth>
+      <superparametric_mapping/>
+    </spherical_earth>
+    <ocean_boundaries>
+      <top_surface_ids>
+        <integer_value shape="1" rank="1">1</integer_value>
+      </top_surface_ids>
+      <bottom_surface_ids>
+        <integer_value shape="1" rank="1">2</integer_value>
+      </bottom_surface_ids>
+      <scalar_field name="DistanceToTop" rank="0">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat/>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </scalar_field>
+      <scalar_field name="DistanceToBottom" rank="0">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat/>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </scalar_field>
+    </ocean_boundaries>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0.0</real_value>
+        <comment>43200: dump results every 12 hours (0.5 days) of simulated time. Given a finish-time of 2592000 s (1 month), 60 dumps are expected. Given a time-step size of 300s, a dump every 144 time steps is expected.</comment>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <checkpointing>
+      <checkpoint_period_in_dumps>
+        <integer_value rank="0">10</integer_value>
+      </checkpoint_period_in_dumps>
+      <checkpoint_at_end/>
+    </checkpointing>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+      <time_units date="seconds since 1987-01-05 00:00.0"/>
+    </current_time>
+    <timestep>
+      <real_value rank="0">3600.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">2592000</real_value>
+      <comment>(2592000)1 month, in seconds. Given a timestep size of 300 s, the simulation should perform 8640 time-steps.</comment>
+    </finish_time>
+    <final_timestep>
+      <integer_value rank="0">2</integer_value>
+    </final_timestep>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+    </nonlinear_iterations>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.81</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X, t):
+
+  a = X[0]
+  b = X[1]
+  c = X[2]
+  
+  x_component = -(a/((a**2+b**2+c**2)**0.5)) 
+  y_component = -(b/((a**2+b**2+c**2)**0.5))
+  z_component = -(c/((a**2+b**2+c**2)**0.5))
+  
+  return [x_component, y_component, z_component]</string_value>
+            </python>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+    <coriolis>
+      <on_sphere>
+        <omega>
+          <real_value rank="0">7.27220522e-5</real_value>
+        </omega>
+      </on_sphere>
+    </coriolis>
+  </physical_parameters>
+  <material_phase name="Fields">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">1.0</real_value>
+          </reference_density>
+          <subtract_out_hydrostatic_level/>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">only first timestep</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="mg">
+            <vertical_lumping/>
+          </preconditioner>
+          <relative_error>
+            <real_value rank="0">1.0e-8</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">5000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <free_surface>
+            <python>
+              <string_value type="code" lines="20" language="python">def val(X,t):
+  from math import asin, atan2, sqrt, pi, exp
+  r=sqrt(sum([x**2 for x in X]))
+  lon=180/pi*atan2(X[1],X[0])
+  lat=180/pi*asin(X[2]/r)
+  d=(lat-52.75)**2+(lon+5.875)**2
+  return exp(-d)</string_value>
+            </python>
+          </free_surface>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <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="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <project_velocity_to_continuous>
+                <mesh name="BaseMesh"/>
+              </project_velocity_to_continuous>
+              <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>
+          <discontinuous_galerkin>
+            <maximum_courant_number_per_subcycle>
+              <real_value rank="0">0.1</real_value>
+            </maximum_courant_number_per_subcycle>
+          </discontinuous_galerkin>
+        </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-8</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</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="FreeSurface">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="free_surface"/>
+        </boundary_conditions>
+        <boundary_conditions name="NoNormalFlow">
+          <surface_ids>
+            <integer_value shape="5" rank="1">2 1009 1010 1011 1013</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_symmetric>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X, t):
+    a = X[0]
+    b = X[1]
+    c = X[2]
+    from math import sqrt, sin, cos, atan2, acos
+    r=sqrt(a**2+b**2+c**2)
+    #A1=2000.0
+    A1=20000.0
+    A2=A1
+    A3=1.0
+    phi=atan2(b,a)
+    theta=acos(c/r)
+    T11=A1*sin(phi)**2+A2*cos(theta)**2*cos(phi)**2+A3*sin(theta)**2*cos(phi)**2
+    T12=-A1*sin(phi)*cos(phi)+A2*cos(theta)**2*sin(phi)*cos(phi)+A3*sin(theta)**2*sin(phi)*cos(phi)
+    T13=-A2*sin(theta)*cos(theta)*cos(phi)+A3*sin(theta)*cos(theta)*cos(phi)
+    T21=T12
+    T22=A1*cos(phi)**2+A2*cos(theta)**2*sin(phi)**2+A3*sin(theta)**2*sin(phi)**2
+    T23=-A2*sin(theta)*cos(theta)*sin(phi)+A3*sin(theta)*cos(theta)*sin(phi)
+    T31=T13
+    T32=T23
+    T33=A2*sin(theta)**2+A3*cos(theta)**2
+    
+    return [[T11, T12, T13],
+           [T21, T22, T23],
+           [T31, T32, T33]]</string_value>
+                </python>
+              </anisotropic_symmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_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="FreeSurface" rank="0">
+      <diagnostic>
+        <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>
+    <scalar_field name="DG_CourantNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="Partition" rank="0">
+      <diagnostic>
+        <algorithm name="element_ownership" material_phase_support="single"/>
+        <mesh name="P0Mesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>
 
=== added directory 'tests/spherical_patch/src'
=== added file 'tests/spherical_patch/src/spherical_segment.geo'
--- tests/spherical_patch/src/spherical_segment.geo	1970-01-01 00:00:00 +0000
+++ tests/spherical_patch/src/spherical_segment.geo	2012-05-16 11:01:22 +0000
@@ -0,0 +1,113 @@
+Printf("Reading function and point definitions...");
+
+Function DrawParallel
+ alpha = parallelSectionStartingY/parallelSectionStartingX;
+ parallelRadius = Sqrt(parallelSectionStartingX^2 + parallelSectionStartingY^2);
+ deltaAlpha = (parallelSectionEndingY/parallelSectionEndingX - parallelSectionStartingY/parallelSectionStartingX)/pointsOnParallel;
+ For t In {1:pointsOnParallel}
+  alpha += deltaAlpha;
+  newParallelPointX = (parallelSectionStartingX/Fabs(parallelSectionStartingX))*parallelRadius/(Sqrt(alpha^2 + 1));
+  newParallelPointY = newParallelPointX*alpha;
+  newPointOnParallel = newp; Point(newPointOnParallel) = {newParallelPointX, newParallelPointY, 0};
+ EndFor
+ BSpline(newParallelID) = {firstPointOnParallel, newPointOnParallel-(pointsOnParallel-1):newPointOnParallel, lastPointOnParallel};
+Return
+
+//The coordinates of various points on the Surface of the Earth are converted to radians 
+//and stereographic coordinatesbelow.
+
+//Point A is located 2 deg 00'W, 48 deg 45'N.
+Point_A_longitude_rad = -(2 + (00/60))*(Pi/180);
+Point_A_latitude_rad = (48 + (45/60))*(Pi/180);
+Point_A_stereographic_x = -Cos(Fabs(Point_A_longitude_rad))*Cos(Fabs(Point_A_latitude_rad))/( 1 + Sin(Fabs(Point_A_latitude_rad)) );
+Point_A_stereographic_y = Cos(Fabs(Point_A_latitude_rad))*Sin(Fabs(Point_A_longitude_rad))/( 1 + Sin(Fabs(Point_A_latitude_rad)) );
+
+//Point B is located at 9 deg 46'W, 48 deg 45'N
+Point_B_longitude_rad = -(9 + (46/60))*(Pi/180);
+Point_B_latitude_rad = (48 + (45/60))*(Pi/180);
+Point_B_stereographic_x = -Cos(Fabs(Point_B_longitude_rad))*Cos(Fabs(Point_B_latitude_rad))/( 1 + Sin(Fabs(Point_B_latitude_rad)) );
+Point_B_stereographic_y = Cos(Fabs(Point_B_latitude_rad))*Sin(Fabs(Point_B_longitude_rad))/( 1 + Sin(Fabs(Point_B_latitude_rad)) );
+
+//Point C is located at 9 deg 45'W, 56 deg 45'N
+Point_C_longitude_rad = -(9 + (45/60))*(Pi/180);
+Point_C_latitude_rad = (56 + (45/60))*(Pi/180);
+Point_C_stereographic_x = -Cos(Fabs(Point_C_longitude_rad))*Cos(Fabs(Point_C_latitude_rad))/( 1 + Sin(Fabs(Point_C_latitude_rad)) );
+Point_C_stereographic_y = Cos(Fabs(Point_C_latitude_rad))*Sin(Fabs(Point_C_longitude_rad))/( 1 + Sin(Fabs(Point_C_latitude_rad)) );
+
+//Point D is located at 2 deg 00'W, 56 deg 45'N
+Point_D_longitude_rad = -(2 + (00/60))*(Pi/180);
+Point_D_latitude_rad = (56 + (45/60))*(Pi/180);
+Point_D_stereographic_x = -Cos(Fabs(Point_D_longitude_rad))*Cos(Fabs(Point_D_latitude_rad))/( 1 + Sin(Fabs(Point_D_latitude_rad)) );
+Point_D_stereographic_y = Cos(Fabs(Point_D_latitude_rad))*Sin(Fabs(Point_D_longitude_rad))/( 1 + Sin(Fabs(Point_D_latitude_rad)) );
+
+//Create point, line, line-loop surface and field ID's
+IP = newp;
+IL = newl;
+ILL = newll;
+IS = news;
+IFI = newf;
+
+//Draw a point on the centre of the earth and one on the North pole, then create a sphere using
+// these points, note that remaining points will be drawin in the stereographic projection plane.
+Point ( IP + 0 ) = {0, 0, 0 };
+Point ( IP + 1 ) = {0, 0,6.37101e+06};
+PolarSphere ( IS + 0 ) = {IP , IP+1};
+
+Printf("Creating Meridian & Parallel segments to close the domain...");
+//So far, only points on the shorelies and curves approximating shorelines have been
+//constructed. we now add points in the sea, to close the domain boundaries
+Point(IP + 60000) = {Point_A_stereographic_x, Point_A_stereographic_y, 0};
+Point(IP + 60001) = {Point_B_stereographic_x, Point_B_stereographic_y, 0};
+Point(IP + 60002) = {Point_C_stereographic_x, Point_C_stereographic_y, 0};
+Point(IP + 60003) = {Point_D_stereographic_x, Point_D_stereographic_y, 0};
+
+//Draw South-most parallel of the Domain.
+pointsOnParallel = 200;               //Assign parameters to variables and then call function DrawParallel,
+parallelSectionStartingX = Point_A_stereographic_x;
+parallelSectionStartingY = Point_A_stereographic_y;
+firstPointOnParallel = IP + 60000;
+parallelSectionEndingX = Point_B_stereographic_x;
+parallelSectionEndingY = Point_B_stereographic_y;
+lastPointOnParallel = IP + 60001;
+newParallelID = IL + 1000;
+Call DrawParallel;
+
+BSpline(IL + 1001) = {IP + 60001, IP + 60002};
+BSpline(IL + 1002) = {IP + 60000, IP + 60003};
+
+////Draw North-most parallel of the Domain.
+pointsOnParallel = 200;            //Assign parameters to variables and then call function DrawParallel,
+parallelSectionStartingX = Point_C_stereographic_x;
+parallelSectionStartingY = Point_C_stereographic_y;
+firstPointOnParallel = IP + 60002;
+parallelSectionEndingX = Point_D_stereographic_x;
+parallelSectionEndingY = Point_D_stereographic_y;
+lastPointOnParallel = IP + 60003;
+newParallelID = IL + 1003;
+Call DrawParallel;
+
+Printf("Declaring surface on sphere to be meshed...");
+Line Loop(1006) = {1003, -1004, -1002, -1001};
+Plane Surface(1007) = {1006};
+Physical Surface(1008) = {1007};
+
+Physical Line(1009) = {1001}; //Physical line corresponding to south-most parallel.
+Physical Line(1010) = {1002}; //Physical line corresponding to meridian segment through point A.
+Physical Line(1011) = {1003}; //Physical line corresponding to meridian segment through point B.
+Physical Line(1013) = {1004}; //Physical line corresponding to north-most parallel.
+
+Printf("Assigning characteristic mesh sizes...");
+Field[ IFI + 1] = MathEval;
+Field[ IFI + 1].F = "50000";
+
+Background Field = IFI + 1;
+
+// Dont extent the elements sizes from the boundary inside the domain
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+
+//Set some options for better png output
+General.Color.Background = {255,255,255};
+General.Color.BackgroundGradient = {255,255,255};
+General.Color.Foreground = Black;
+Mesh.Color.Lines = {0,0,0};
+
 
=== modified file 'tools/Checkmesh.F90'
--- tools/Checkmesh.F90	2012-04-03 12:24:30 +0000
+++ tools/Checkmesh.F90	2012-05-16 11:01:22 +0000
@@ -15,7 +15,7 @@
   use linked_lists
   use meshdiagnostics
   use metric_tools
-  use mesh_files
+  use read_triangle
   use supermesh_construction
   use tetrahedron_intersection_module
   use iso_c_binding
@@ -38,7 +38,7 @@
   rformat = real_format()    
 
   print "(a)", "Reading in mesh mesh with base name " // trim(filename)
-  positions = read_mesh_files(filename, quad_degree = 4)    
+  positions = read_triangle_files(filename, quad_degree = 4)    
   if(isparallel()) call read_halos(filename, positions)
   print "(a)", "Read successful"
 
 
=== modified file 'tools/Fladapt.F90'
--- tools/Fladapt.F90	2012-04-16 12:54:50 +0000
+++ tools/Fladapt.F90	2012-05-16 11:01:22 +0000
@@ -76,6 +76,7 @@
   type(vector_field) :: new_mesh_field
   type(vector_field), pointer :: new_mesh_field_ptr, old_mesh_field
   type(tensor_field) :: metric, t_edge_lengths
+  character(len=FIELD_NAME_LEN) :: mesh_format
 
   ! now turn into proper fortran strings (is there an easier way to do this?)
   do i=1, input_basename_len
@@ -144,11 +145,10 @@
   else
     call adapt_mesh(old_mesh_field, metric, new_mesh_field)
   end if
-  old_mesh_field => null()
-  old_mesh => null()
   
+  call get_option(trim(old_mesh_field%mesh%option_path)//"/from_file/format/name", mesh_format)
   ! Write the output mesh
-  call write_mesh_files(output_basename, new_mesh_field)
+  call write_mesh_files(output_basename, mesh_format, new_mesh_field)
   
   ! Deallocate
   do i = 1, size(states)
@@ -183,7 +183,7 @@
       FLAbort("External mesh field " // trim(mesh_field_name) // " not found in the system state")
     end if
       
-    mesh_field => extract_vector_field(state, mesh_field_name)    
+    mesh_field => extract_vector_field(state, mesh_field_name)
     
   end subroutine find_mesh_field_to_adapt
   
 
=== modified file 'tools/periodise.F90'
--- tools/periodise.F90	2011-02-28 14:43:05 +0000
+++ tools/periodise.F90	2012-05-16 11:01:22 +0000
@@ -20,7 +20,7 @@
   integer :: status
   type(state_type), dimension(:), pointer :: states
   integer :: ierr
-  character(len=FIELD_NAME_LEN) :: external_name, periodic_name
+  character(len=FIELD_NAME_LEN) :: external_name, periodic_name, mesh_format
   type(mesh_type), pointer :: periodic_mesh, external_mesh
   type(vector_field) :: periodic_positions, external_positions
   integer :: stat, i, nstates
@@ -83,7 +83,8 @@
 
   ! Dump out the periodic mesh to disk:
   new_external_filename = trim(external_filename) // '_periodic'
-  call write_mesh_files(new_external_filename, periodic_positions)
+  call get_option(trim(external_mesh%option_path)//"/from_file/format/name", mesh_format)
+  call write_mesh_files(new_external_filename, mesh_format, periodic_positions)
 
   ! OK! Now we need to do some setting of options.
   call manipulate_options(external_mesh, trim(external_mesh%option_path), periodic_mesh, trim(periodic_mesh%option_path), new_external_filename)
 
=== modified file 'tools/project_to_continuous.F90'
--- tools/project_to_continuous.F90	2011-09-09 13:24:20 +0000
+++ tools/project_to_continuous.F90	2012-05-16 11:01:22 +0000
@@ -35,7 +35,7 @@
   use state_module
   use elements
   use fields
-  use mesh_files
+  use read_triangle
   use vtk_interfaces
   use sparse_tools
   use fefields
@@ -72,7 +72,7 @@
 
   call vtk_read_state(vtuname, dg_state, quad_degree=6)
   
-  cg_coordinate= read_mesh_files(meshname, quad_degree=6)
+  cg_coordinate= read_triangle_files(meshname, quad_degree=6)
   cg_mesh=cg_coordinate%mesh
 
   call allocate(lumped_mass, cg_mesh, "LumpedMass")
 
=== modified file 'tools/test_laplacian.F90'
--- tools/test_laplacian.F90	2010-07-02 08:26:03 +0000
+++ tools/test_laplacian.F90	2012-05-16 11:01:22 +0000
@@ -9,7 +9,7 @@
   ! 1x1x1 cube with the boundary conditions 1 and 2 applied on the
   ! sides of the first coordinates direction (e.g. test_laplacian.poly
   ! - create a mesh with 'triangle -a0.01 -e test_laplacian.poly' )
-  use mesh_files
+  use read_triangle
   use fields
   use FEtools
   use elements
@@ -70,7 +70,7 @@
 
   call read_command_line(filename, degree, quad_degree)
 
-  positions=read_mesh_files(filename, quad_degree=quad_degree)
+  positions=read_triangle_files(filename, quad_degree=quad_degree)
 
   call insert(state, positions, "Coordinate")