Merge lp:~fluidity-core/fluidity/prescribed-regions-adaptivity into lp:fluidity

[Stephan Kramer]

Improved support for internal boundaries, in particular 2d adaptivity.

This branch is a bit of a misnomer really, there is nothing specifically about prescribed regions in this merge. After this merge, you can:

* read in a gmsh or triangle mesh with an internal boundary directly.
Previously a hack was needed: by running gmsh2triangle with an extra -i option, the internal facets would be doubled (as needed for Fluidity's data structures) and an extra column in the .edge/.face file would be written out that indicates which element each facet belongs to, so we can distinguish between the two coincident facets. Having both facets in the input mesh also enables the two facets to have a different surface id, thus allowing different boundary conditions on either side of the internal boundary. This functionality is also used in checkpointing of periodic meshes (in particular in periodise) with an aliased and physical surface id on either side of the periodic boundary, which is an internal boundary in the periodic mesh.
For most other practical purposes however having a different boundary condition on either side is not really possible due to restrictions of the assembly code (and in the case of cg having only a single degree of freedom). Moreover a mesh prepared in gmsh will only allow you to to assign one surface id, and will only write out a single facet.
Thus, to cut a long story short, this merge allows you to read in a gmsh file with single facets along the internal boundary directly into Fluidity. The facets will be doubled inside add_faces() and will of course be assigned the same surface id.
* use 2d adaptivity with an internal boundary
Here, only the single valued internal boundary condition is supported. This is because the 2d adaptivity library only deals with single faceted internal boundaries. The internal boundary is given as a constraint to the library, and returned (in the new numbering) after the adapt. A call to add_faces() then takes care of the doubling of the internal facets.

To deal with these changes the following interface changes were introduced:
* rename has_internal_boundaries() to has_discontinuous_internal_boundaries()
This is because this routine now only applies to the case where the internal boundary does have multi-valued surface ids and is read and written with an element_owner column in the mesh file. At the moment this is only used for periodic meshes.
* add unique_surface_element_count()
When the internal facets are duplicated inside add_faces(), the new ones are added at the end of the list of the original, unique facets provided in the mesh. surface_element_count() will include both copies and should be used in most places in the assembly code dealing with boundary conditions. unique_surface_element_count() only counts the original single facets and is only used in some special cases, such as providing the facets to the adaptivity library, or when writing out the .edge/.face file.

This functionality is tested in a new test, internal_boundary_adaptivity.

NOTE that the actual applying of internal boundaries to the equations is not very well tested. This depends very much on which discretisation is used and which type of bcs. It should therefore be carefully evaluated for the case you want to apply it to. You will probably have most luck with strong bcs, with weak bcs your mileage may vary.

[Rhodri Davies (rhodri-davies)]

I have run a suite of 2D simulations, in serial and parallel, to test this branch today. They all worked beautifully. A glance at the code suggests that everything is nicely done, although I suspect there may be somebody better to look at certain aspects than myself. I note that in adapt integration you have the warning:

    if (surface_element_count(input_positions)/=unique_surface_element_count(input_positions%mesh)) then
      ewrite(0,*) "WARNING: It appears you have an internal boundary and you're trying to use 3D adaptivity."
      ewrite(0,*) "This combination has not been tested yet and may not work at all."
      ewrite(0,*) "Please report your experiences to the Fluidity mailing list."
    end if

If you're happy to wait a few days, I can test this with one of Giuseppe's 3D cases and feedback. However, if you're keen to get this merged I can easily do that post merge. Your choice entirely.

Rhod

[Stephan Kramer (s-kramer)]

Thanks for testing this, Rhodri!

That warning exists not because "in theory it should work but I haven't had time to test it". It's simply not implemented, the stuff I added for 2d adaptivity is in a completely different code path. So there's three possible outcomes:
1) It crashes immediately
2) It goes through an adapt but the internal boundary dissapears (this is what would happen in 2D before)
3) It somehow magically works (quite unlikely)

Now if there's demand I'm quite happy at implementing it for 3d as well (this should be a lot easier now with the new add_faces() functionality), but I don't think this needs to go in this merge.

[Jon Hill (jon-hill)]

 - Can you reword the warning from "not tested" to "not implemented". This is more accurate
 - "owner ship" is one word (Fields_Allocates)

Otherwise looks great. I'll leave Sam to do the approve.

[lp:~fluidity-core/fluidity/prescribed-regions-adaptivity]

4229. By Stephan Kramer on 2013-10-16

Merge from trunk.

4230. By Stephan Kramer on 2013-10-16

Changes to error messages, as suggested by jhill.

[Samuel Parkinson (s-parkinson11)]

Hey Stephan,

Generally all top notch. Just a few comments:
1) in Mba2D_integration.F90 there is a comment saying "not sure this will work with lock_faces". What does this mean?
2) lots of testing which is good, but errors get displayed when running the new test saying "The following boundary ids were specified, but they don't appear in the surface mesh:".
3) You mention that only certain discretisations are supported. Should there be a manual entry listing the caveats under which this is tested and is known to work.

Other than that it's all good.

[Stephan Kramer (s-kramer)]

Hey Sam,

Thanks for the review.
Ad 1): had to dig deep in my memory. It seems lock_faces is only used with periodic adaptivity. I think in general internal boundaries would need testing with that anyway. What I specifically had in mind with that comment I can't really work out, so I just removed it.
Ad 2): These are spurious warnings actually. This is a separate problem with prescribed regions. It calls populate state on the substate (the state that is based on the mesh with the prescribed regions removed) and then finds that some of the boundaries do no exist. I have opened a bug report (https://bugs.launchpad.net/fluidity/+bug/1240921) as that should be fixed.
Ad 3): Good point. Will add some stuff to the manual with a lot of caveats.

[lp:~fluidity-core/fluidity/prescribed-regions-adaptivity]

4231. By Stephan Kramer on 2013-10-17

Remove unclear comment.

4232. By Stephan Kramer on 2013-10-17

Say something about internal bcs in the manual, as suggested by sparkinson.

[Samuel Parkinson (s-parkinson11)]

review: Approve

=== modified file 'assemble/Adapt_Integration.F90'
--- assemble/Adapt_Integration.F90	2011-04-22 17:28:48 +0000
+++ assemble/Adapt_Integration.F90	2013-10-17 12:55:48 +0000
@@ -357,6 +357,15 @@
     if(nselm > 0) then
       call getsndgln(input_positions%mesh, snlist)
     end if
+
+    if (surface_element_count(input_positions)/=unique_surface_element_count(input_positions%mesh)) then
+      ewrite(0,*) "It appears you have an internal boundary and you're trying to use 3D adaptivity."
+      ewrite(0,*) "This combination has not been implemented yet."
+      ! You could try to see if it somehow does work, by simply removing this FLExit()
+      ! (make sure to check you still have the right internal boundary ids after the adapt)
+      ! Feel free to discuss on the fluidity mailing list.
+      FLExit("Cannot have internal boundaries with 3D adaptivity")
+    end if
     
     ! Surface IDs
     allocate(surfid(nselm))

=== modified file 'assemble/Mba2d_Integration.F90'
--- assemble/Mba2d_Integration.F90	2010-11-25 12:09:32 +0000
+++ assemble/Mba2d_Integration.F90	2013-10-17 12:55:48 +0000
@@ -47,7 +47,7 @@
 
     type(mesh_type), pointer :: xmesh
 
-    integer :: nonods, mxnods, orig_stotel, stotel, mxface, totele, maxele
+    integer :: nonods, mxnods, orig_stotel, stotel, mxface, totele, maxele, stotel_external
     real, dimension(:, :), allocatable :: pos
     integer, dimension(:, :), allocatable :: ipf
     integer, dimension(:, :), allocatable :: ipe
@@ -57,7 +57,7 @@
     integer, dimension(:), allocatable :: iFnc
     integer, dimension(:), allocatable :: lbE
     real, dimension(:, :), allocatable :: tmp_metric
-    integer :: i, j, k, partition_surface_id, face
+    integer :: i, j, k, partition_surface_id, face, face2
     real :: quality, rQuality
     integer :: iPrint, ierr, maxWr, maxWi
     real, dimension(:), allocatable :: rW
@@ -144,6 +144,7 @@
     ipf = 0
     partition_surface_id = maxval(surface_ids) + 1
     stotel = 0
+    stotel_external = 0
 
     call allocate(input_face_numbering_to_mba2d_numbering)
     do i=1,totele
@@ -156,8 +157,9 @@
           call insert(input_face_numbering_to_mba2d_numbering, face, stotel)
           ipf(1:2, stotel) = face_global_nodes(xmesh, face)
           ipf(3, stotel) = 0
-          if (face <= orig_stotel) then ! if face is genuinely external, i.e. on the domain exterior
+          if (face <= orig_stotel) then ! if facet is genuinely external, i.e. on the domain exterior
             ipf(4, stotel) = surface_ids(face)
+            stotel_external = stotel_external + 1
           else
             ipf(4, stotel) = partition_surface_id
           end if
@@ -165,6 +167,28 @@
       end do
     end do
 
+    if (stotel_external<orig_stotel) then
+      ! not all facets in the surface mesh are external apparently
+      do face=1, orig_stotel
+        face2 = face_neigh(xmesh, face)
+        if (face>face2) then
+          ! if face==face2 it's an external facet that's dealt with already
+          ! we check face>face2 to ensure we only copy one of the two coinciding facets
+          stotel = stotel +1
+          call insert(input_face_numbering_to_mba2d_numbering, face, stotel)
+          ipf(1:2, stotel) = face_global_nodes(xmesh, face)
+          ipf(3, stotel) = 0
+          ipf(4, stotel) = surface_ids(face)
+          if (surface_ids(face)/=surface_ids(face2)) then
+            ! note that we're here checking for the combined surface+coplanar id to be the same
+            ! but the coplanar id of the 2 facets should always be the same
+            FLExit("Adaptivity with internal boundaries only works if the surface id is single valued")
+          end if
+        end if                
+      end do
+
+    end if
+
     if (.not. present(lock_faces)) then
       nfv = 0
       allocate(ifv(nfv))
@@ -343,22 +367,9 @@
     new_mesh%option_path = xmesh%option_path
     
     if (.not. isparallel()) then
-      allocate(boundary_ids(stotel))
-      allocate(coplanar_ids(stotel))
-      call deinterleave_surface_ids(ipf(4, 1:stotel), max_coplanar_id, boundary_ids, coplanar_ids)  
-      
-      allocate(new_sndgln(1:2,1:stotel))
-      new_sndgln=ipf(1:2,1:stotel)
-      call add_faces(new_mesh, sndgln=reshape(new_sndgln, (/ 2*stotel /) ), &
-        boundary_ids=boundary_ids)
-      deallocate(boundary_ids)
-      deallocate(new_sndgln)
-      
-      if(associated(xmesh%faces%coplanar_ids)) then
-        allocate(new_mesh%faces%coplanar_ids(stotel))
-        new_mesh%faces%coplanar_ids = coplanar_ids
-      end if
-      deallocate(coplanar_ids)
+      allocate(new_sndgln(1:2,1:stotel), mba_boundary_ids(1:stotel))
+      new_sndgln = ipf(1:2,1:stotel)
+      mba_boundary_ids =ipf(4,1:stotel)
     else
       ! In parallel, we need to filter out the surface elements with colour partition_surface_id, because
       ! they are not real external faces
@@ -372,35 +383,39 @@
         end if
       end do
 
+      ! number of surface elements without inter-partition surface elements
       stotel = j - 1
-      allocate(mba_boundary_ids(stotel))
-
-      do i=1,stotel
+      allocate(mba_boundary_ids(1:stotel), new_sndgln(1:2,1:stotel))
+      do i=1, stotel
         mba_boundary_ids(i) = ipf(4, fetch(physical_surface_ids, i))
-      end do
-
-      allocate(boundary_ids(stotel))
-      allocate(coplanar_ids(stotel))
-      call deinterleave_surface_ids(mba_boundary_ids, max_coplanar_id, boundary_ids, coplanar_ids)  
-      deallocate(mba_boundary_ids)
-      
-      allocate(new_sndgln(1:2,1:stotel))
-      do i=1,stotel
-        new_sndgln(1:2, i) = ipf(1:2, fetch(physical_surface_ids, i))
-      end do
-      call add_faces(new_mesh, sndgln=reshape(new_sndgln, (/ 2*stotel /) ), &
-        boundary_ids=boundary_ids)
-      deallocate(boundary_ids)
-      deallocate(new_sndgln)
-      
-      if(associated(xmesh%faces%coplanar_ids)) then
-        allocate(new_mesh%faces%coplanar_ids(stotel))
-        new_mesh%faces%coplanar_ids = coplanar_ids
-      end if
-      deallocate(coplanar_ids)
-
+        new_sndgln(1:2, i) = ipf(1:2, fetch(physical_surface_ids, i))
+      end do
+      
+      do i=1, stotel
+        new_sndgln(1:2, i) = ipf(1:2, fetch(physical_surface_ids, i))
+      end do
       call deallocate(physical_surface_ids)
     end if
+
+    ! add_faces might create extra (internal) surface elements, so we 
+    ! use the combined boundary+coplanar ids first
+    call add_faces(new_mesh, sndgln=reshape(new_sndgln, (/ 2*stotel /) ), &
+      boundary_ids=mba_boundary_ids)
+    deallocate(mba_boundary_ids)
+    deallocate(new_sndgln)
+
+    ! and only deinterleave now we know the total number of elements in the surface mesh
+    stotel = surface_element_count(new_mesh)
+    allocate(boundary_ids(1:stotel), coplanar_ids(1:stotel))
+    call deinterleave_surface_ids(new_mesh%faces%boundary_ids, max_coplanar_id, boundary_ids, coplanar_ids)
+
+    new_mesh%faces%boundary_ids = boundary_ids
+
+    if(associated(xmesh%faces%coplanar_ids)) then
+      allocate(new_mesh%faces%coplanar_ids(1:stotel))
+      new_mesh%faces%coplanar_ids = coplanar_ids
+    end if
+    deallocate(boundary_ids, coplanar_ids)
     
     if(have_option("/mesh_adaptivity/hr_adaptivity/preserve_mesh_regions")&
                               .or.present_and_true(force_preserve_regions)) then

=== modified file 'femtools/Adjacency_Lists.F90'
--- femtools/Adjacency_Lists.F90	2010-11-09 14:08:19 +0000
+++ femtools/Adjacency_Lists.F90	2013-10-17 12:55:48 +0000
@@ -908,7 +908,7 @@
            if(.not. any(row_idx(j - 1)%ptr == candidate_ele)) cycle ele_loop
          end do
          
-#ifndef DDEBUG
+#ifdef DDEBUG
          if(adj_ele > 0) then
            ! We've found more than one adjacent element. We're in trouble.
            adj_ele = -1

=== modified file 'femtools/Fields_Allocates.F90'
--- femtools/Fields_Allocates.F90	2013-03-05 10:21:07 +0000
+++ femtools/Fields_Allocates.F90	2013-10-17 12:55:48 +0000
@@ -154,7 +154,7 @@
     integer, intent(in) :: nodes, elements
     type(element_type), target, intent(in) :: shape
     character(len=*), intent(in), optional :: name
-    integer :: i, j
+    integer :: i
     
     mesh%nodes=nodes
 
@@ -909,7 +909,7 @@
     
     integer, dimension(:), allocatable :: ndglno
     real, dimension(:), pointer :: val
-    integer :: i, j, input_nodes, n_faces
+    integer :: i, input_nodes, n_faces
 
     if (present(continuity)) then
        mesh%continuity=continuity
@@ -1126,9 +1126,15 @@
     !! be consistent with that of the periodic face local node numbering.
     type(integer_hash_table), intent(in), optional :: periodic_face_map
     !! This list gives the element owning each face in sndgln. This needs
-    !! to be provided when sndgln contains internal faces, as we need to
-    !! decide which of the two adjacent elements the face belongs to
-    !! (used in reading in periodic meshes which include the periodic faces in the surface mesh)
+    !! to be provided when sndgln contains internal faces and both copies
+    !! are included, as we need to decide which of the two adjacent elements 
+    !! the facet belongs to (used in reading in periodic meshes which include 
+    !! the periodic facets in the surface mesh with a different surface id on 
+    !! either side of the periodic boundary). If no element_owner information
+    !! is provided, internal facets in sndgln are assumed to only appear once
+    !! and a copy will be made, its boundary id will be copied as well. This 
+    !! means afterwards the surface_element_count() will be higher than the 
+    !! number of facets provided in sndgln.
     integer, dimension(:), intent(in), optional :: element_owner
     !! See comments above sndgln
     logical, intent(in), optional :: incomplete_surface_mesh
@@ -1211,8 +1217,8 @@
          ! works as long as we're not discontinuous
          assert( mesh%continuity>=0 )
          
-         ! the periodic faces will be internal faces in the output periodic mesh
-         mesh%faces%has_internal_boundaries = .true.         
+         ! the periodic faces will be discontinuous internal faces in the output periodic mesh
+         mesh%faces%has_discontinuous_internal_boundaries = .true.         
          
       else if (model%periodic .and. .not. mesh%periodic) then
         
@@ -1221,7 +1227,7 @@
             mesh, model, lperiodic_face_map, stat=stat)
             
          ! the subroutine above only works if the removing of periodic bcs has removed all internal boundaries
-         mesh%faces%has_internal_boundaries = .true.
+         mesh%faces%has_discontinuous_internal_boundaries = .false.
                      
       else
          ! Transfer the faces from model to mesh
@@ -1229,7 +1235,7 @@
          call incref(mesh%faces%face_list)
          
          ! have internal faces if the model does
-         mesh%faces%has_internal_boundaries = has_internal_boundaries(model)
+         mesh%faces%has_discontinuous_internal_boundaries = has_discontinuous_internal_boundaries(model)
       end if
         
       ! face_element_list is a pure copy of that of the model
@@ -1249,6 +1255,7 @@
            size(mesh%faces%boundary_ids), &
            trim(mesh%name)//" boundary_ids")
 #endif
+      mesh%faces%unique_surface_element_count = model%faces%unique_surface_element_count
        
       ! same for coplanar ids (if existent)
       if (associated(model%faces%coplanar_ids)) then
@@ -1382,15 +1389,16 @@
     integer, dimension(:), intent(in), optional :: element_owner
     logical, intent(in), optional :: incomplete_surface_mesh
 
+    type(integer_hash_table):: internal_facet_map
     type(element_type), pointer :: mesh_shape
     type(element_type) :: face_shape
-    logical:: internal_face, surface_elements_added
+    logical:: surface_elements_added
     ! listen very carefully, I shall say this only once:
     logical, save :: warning_given=.false. 
     integer, dimension(:), pointer :: faces, neigh, snodes
     integer, dimension(2) :: common_elements
-    integer :: nloc, snloc, stotel
-    integer :: bdry_count, ele, sele, j, no_found
+    integer :: snloc, stotel
+    integer :: bdry_count, ele, sele, sele2, neighbour_ele, j, no_found
     integer :: no_faces
     type(csr_sparsity) :: sparsity
     type(csr_sparsity) :: nelist
@@ -1398,7 +1406,6 @@
     assert(continuity(mesh)>=0)
 
     mesh_shape=>ele_shape(mesh, 1)
-    nloc=mesh_shape%loc
 
     ! Calculate the node-to-element list.
     ! Calculate the element adjacency list.
@@ -1415,12 +1422,16 @@
          trim(mesh%name)//" face_element_list")
 #endif
 
-    mesh%faces%has_internal_boundaries = present(element_owner)
+    mesh%faces%has_discontinuous_internal_boundaries = present(element_owner)
+
+    call allocate(internal_facet_map)
     
     snloc=face_vertices(mesh_shape)
     if (present(sndgln)) then
       
-      stotel=size(sndgln)/snloc
+      stotel=size(sndgln)/snloc ! number of provided surface elements
+      ! we might add some more when duplicating internal facets
+      bdry_count=stotel
       
       do sele=1, stotel
         
@@ -1432,63 +1443,67 @@
         if (no_found==1) then
           ! we have found the adjacent element
           ele=common_elements(1)
-          internal_face=.false.
           if (present(element_owner)) then
-            assert(element_owner(sele)==ele)
+            if (element_owner(sele)/=ele) then
+              ewrite(0,*) "Surface element: ", sele
+              ewrite(0,*) "Provided element owner: ", element_owner(sele)
+              ewrite(0,*) "Found adjacent element: ", ele
+              FLExit("Provided element ownership information is incorrect")
+            end if
           end if
+          call register_external_surface_element(mesh, sele, ele, snodes)
         else if (no_found==2 .and. present(element_owner)) then
-          ele=element_owner(sele)
-          if (all(common_elements/=ele)) then
-            FLAbort("Face not adjacent to specified element_owner")
-          end if
-          internal_face=.true.
-        else
-          FLAbort("Surface element in sndgln is not on the surface.")
-        end if
-        
-        mesh%faces%face_element_list(sele)=ele
-        
-        ! now we have to fill in face_list:
-        neigh=>row_m_ptr(mesh%faces%face_list, ele)
-        faces=>row_ival_ptr(mesh%faces%face_list, ele)
-
-        ! find the matching boundary of this element
-        do j=1, mesh%shape%numbering%boundaries
-          if (neigh(j)<=0 .or. internal_face) then
-            if (SetContains(snodes, mesh%ndglno( (ele-1)*nloc+ &
-                                  boundary_numbering(mesh%shape%numbering, j) &
-                                            ))) exit
-          end if
-        end do
-          
-        if (j>mesh%shape%numbering%boundaries) then
-          ! not found a matching boundary, something's wrong
-          FLAbort("Something wrong with the mesh, sndgln, or mesh%nelist")
-          
-        else if (neigh(j)/=0 .and. .not. internal_face) then
-          ! this surface element is already registered
-          ! so apparently there's a duplicate element in the surface mesh
-          ewrite(0,*) 'Surface element:', faces(j),' and ',sele
-          ewrite(0,*) 'Both define the surface element:', snodes
-          FLAbort("Duplicate element in the surface mesh")
-        end if
-        
-        ! register the surface element in face_list
-        faces(j)=sele
-        if (.not. internal_face) then
-          neigh(j)=-j ! negative number indicates exterior boundary
-        else
-          ! if all went well neigh(j) should have the right value:
-          assert( neigh(j)==minval(common_elements, mask=common_elements/=ele) )
+          ! internal facet with element ownership infomation provided
+          ! so we assume both of the coinciding internal facets are
+          ! present in the provided surface mesh and register only
+          ! one of them here
+          ele = element_owner(sele)
+          if (ele==common_elements(1)) then
+            neighbour_ele = common_elements(2)
+          else if (ele==common_elements(2)) then
+            neighbour_ele = common_elements(1)
+          else
+            ewrite(0,*) "Surface element: ", sele
+            ewrite(0,*) "Provided element owner: ", ele
+            ewrite(0,*) "Found adjacent elements: ", common_elements
+            FLExit("Provided element owner ship information is incorrect")
+          end if
+          call register_internal_surface_element(mesh, sele, ele, neighbour_ele)
+        else if (no_found==2) then
+          ! internal facet but not element ownership information:
+          ! we assume this facet only occurs once and we register both
+          ! copies at once
+
+          ! first one using the current surface element number:
+          call register_internal_surface_element(mesh, sele, common_elements(1), common_elements(2))
+
+          ! for the second one we create a new facet number at the end of the provided number surface
+          ! elements:
+          bdry_count = bdry_count+1
+          call register_internal_surface_element(mesh, bdry_count, common_elements(2), common_elements(1))
+          ! store this pair so we can later copy the boundary id of the first (sele) to the second one (bdry_count)
+          call insert(internal_facet_map, sele, bdry_count)
+
+        else if (no_found==0) then
+          ewrite(0,*) "Current surface element: ", sele
+          ewrite(0,*) "With nodes: ", snodes
+          FLExit("Surface element does not exist in the mesh")
+        else
+          ! no_found>2 apparently
+          ! this might happen when calling add_faces on meshes with non-trivial toplogy such
+          ! as calling add_faces on the surface_mesh - we can't really deal with that
+          ewrite(0,*) "Current surface element: ", sele
+          ewrite(0,*) "With nodes: ", snodes
+          FLExit("Surface element (facet) adjacent to more than two elements!")
         end if
         
       end do
         
-      bdry_count=stotel
       
     else
     
-      bdry_count=0
+      stotel=0 ! number of facets provided in sndgln
+      bdry_count=0 ! number of facets including the doubling of interior facets
       
     end if
             
@@ -1517,12 +1532,22 @@
          end do
          
       end do
-      if (surface_elements_added .and. .not. warning_given) then
+
+      if (surface_elements_added .and. key_count(internal_facet_map)>0) then
+        ewrite(0,*) "It appears this mesh has internal boundaries."
+        ewrite(0,*) "In this case all external boundaries need to be marked with a surface id."
+        FLExit("Incomplete surface mesh")
+      else if (surface_elements_added .and. .not. warning_given) then
         ewrite(0,*) "WARNING: an incomplete surface mesh has been provided."
         ewrite(0,*) "This will not work in parallel."
         ewrite(0,*) "All parts of the domain boundary need to be marked with a (physical) surface id."
         warning_given=.true.
       end if
+
+      if (surface_elements_added) then
+        stotel = bdry_count
+      end if
+
     end if
       
     ! the size of this array will be the way to store the n/o
@@ -1533,15 +1558,24 @@
          trim(mesh%name)//" boundary_ids")
 #endif
 
+    mesh%faces%unique_surface_element_count = stotel
+    ewrite(2,*) "Number of surface elements: ", bdry_count
+    ewrite(2,*) "Number of unique surface elements: ", stotel
+
     mesh%faces%boundary_ids=0
     ! copy in supplied boundary ids
     if (present(boundary_ids)) then
       if (.not. present(sndgln)) then
         FLAbort("Boundary ids can only be supplied to add_faces with associated surface mesh")
-      else if (stotel/=size(boundary_ids)) then
+      else if (size(boundary_ids) /= size(sndgln)/snloc) then
         FLAbort("Must supply boundary_ids array for the same number of elements as the surface mesh sndgln")
       end if
-      mesh%faces%boundary_ids(1:stotel)=boundary_ids
+      mesh%faces%boundary_ids(1:size(boundary_ids))=boundary_ids
+
+      do j=1, key_count(internal_facet_map)
+        call fetch_pair(internal_facet_map, j, sele, sele2)
+        mesh%faces%boundary_ids(sele2) = boundary_ids(sele)
+      end do
     end if
 
     ! register the rest of the boundaries (the interior ones):
@@ -1572,9 +1606,86 @@
     assert(bdry_count==size(mesh%faces%face_list%sparsity%colm))
     assert(.not.any(mesh%faces%face_list%sparsity%colm==0))
     assert(.not.any(mesh%faces%face_list%ival==0))
+
+    call deallocate(internal_facet_map)
     
   end subroutine add_faces_face_list
+
+  subroutine register_internal_surface_element(mesh, sele, ele, neighbour_ele)
+    type(mesh_type), intent(inout):: mesh
+    integer, intent(in):: sele, ele, neighbour_ele
+
+    integer, dimension(:), pointer:: neigh, faces
+    integer:: j
+
+    mesh%faces%face_element_list(sele)=ele
+
+    ! neigh should contain correct neighbours already for internal facets:
+    neigh=>row_m_ptr(mesh%faces%face_list, ele)
+    faces=>row_ival_ptr(mesh%faces%face_list, ele)
+
+    ! find the corresponding boundary of ele
+    ! by searching for neighbour_ele in neigh
+    do j=1, mesh%shape%numbering%boundaries
+      if (neigh(j)==neighbour_ele) exit
+    end do
+
+    if (j>mesh%shape%numbering%boundaries) then
+      ! not found a matching boundary, something's wrong
+      FLAbort("Something wrong with the mesh, sndgln, or mesh%nelist")
+    end if
+
+    ! register the surface element in face_list
+    faces(j)=sele
+
+  end subroutine register_internal_surface_element
+
+  subroutine register_external_surface_element(mesh, sele, ele, snodes)
+    type(mesh_type), intent(inout):: mesh
+    integer, intent(in):: sele, ele
+    integer, dimension(:), intent(in):: snodes
+
+    integer, dimension(:), pointer:: neigh, faces, nodes
+    integer:: j, nloc
+
+    nloc = ele_loc(mesh, ele)
     
+    mesh%faces%face_element_list(sele)=ele
+
+    ! for external facets the corresponding entry in neigh
+    ! is still zero, a negative value will be filled in
+    ! when it's registered
+    neigh=>row_m_ptr(mesh%faces%face_list, ele)
+    faces=>row_ival_ptr(mesh%faces%face_list, ele)
+    nodes => ele_nodes(mesh, ele)
+
+    ! find the matching boundary of ele
+    do j=1, mesh%shape%numbering%boundaries
+      ! also check negative entries to check for duplicate registrations
+      if (neigh(j)<=0) then
+        if (SetContains(snodes, nodes(boundary_numbering(mesh%shape%numbering, j)))) exit
+      end if
+    end do
+
+    if (j>mesh%shape%numbering%boundaries) then
+      ! not found a matching boundary, something's wrong
+      FLAbort("Something wrong with the mesh, sndgln, or mesh%nelist")
+    end if
+
+    if (neigh(j)/=0) then
+      ! this surface element is already registered
+      ! so apparently there's a duplicate element in the surface mesh
+      ewrite(0,*) 'Surface element:', faces(j),' and ',sele
+      ewrite(0,*) 'Both define the surface element:', snodes
+      FLAbort("Duplicate element in the surface mesh")
+    end if
+
+    ! register the surface element in face_list
+    faces(j)=sele
+    neigh(j)=-j ! negative number indicates exterior boundary
+
+  end subroutine register_external_surface_element
+
   subroutine add_faces_face_list_periodic_from_non_periodic_model( &
      mesh, model, periodic_face_map)
      ! computes the face_list of a periodic mesh by copying it from

=== modified file 'femtools/Fields_Base.F90'
--- femtools/Fields_Base.F90	2013-07-31 13:45:29 +0000
+++ femtools/Fields_Base.F90	2013-10-17 12:55:48 +0000
@@ -223,6 +223,10 @@
           & surface_element_count_tensor, surface_element_count_mesh
   end interface
 
+  interface unique_surface_element_count
+    module procedure unique_surface_element_count_mesh
+  end interface
+
   interface surface_element_id
     module procedure surface_element_id_scalar, surface_element_id_vector, &
       surface_element_id_mesh
@@ -262,9 +266,9 @@
           & mesh_periodic_tensor
   end interface
     
-  interface has_internal_boundaries
-     module procedure mesh_has_internal_boundaries
-  end interface has_internal_boundaries
+  interface has_discontinuous_internal_boundaries
+     module procedure mesh_has_discontinuous_internal_boundaries
+  end interface has_discontinuous_internal_boundaries
 
   interface extract_scalar_field ! extract_scalar_field is already used in State.F90
      module procedure extract_scalar_field_from_vector_field, extract_scalar_field_from_tensor_field
@@ -412,18 +416,21 @@
   
   end function mesh_periodic_tensor
 
-  pure function mesh_has_internal_boundaries(mesh)
-    !!< Return whether the mesh has internal boundaries
-    logical :: mesh_has_internal_boundaries
+  pure function mesh_has_discontinuous_internal_boundaries(mesh)
+    !!< Return whether the mesh has discontinuous boundaries
+    !!< These are internal boundaries where the surface id are 
+    !!< allowed to be discontinuous (the pair of adjacent 
+    !!< internal facets can have two different ids).
+    logical :: mesh_has_discontinuous_internal_boundaries
     type(mesh_type), intent(in) :: mesh
     
     if (associated(mesh%faces)) then
-      mesh_has_internal_boundaries = mesh%faces%has_internal_boundaries
+      mesh_has_discontinuous_internal_boundaries = mesh%faces%has_discontinuous_internal_boundaries
     else
-      mesh_has_internal_boundaries = .false.
+      mesh_has_discontinuous_internal_boundaries = .false.
     end if
   
-  end function mesh_has_internal_boundaries
+  end function mesh_has_discontinuous_internal_boundaries
 
   pure function node_count_mesh(mesh) result (node_count)
     ! Return the number of nodes in a mesh.
@@ -575,6 +582,21 @@
 
   end function surface_element_count_mesh  
   
+  pure function unique_surface_element_count_mesh(mesh) result (element_count)
+    ! Return the number of unique surface elements of a mesh. For internal
+    ! facets that are part of the surface mesh, each pair of adjacent facets
+    ! is only counted once.
+    integer :: element_count
+    type(mesh_type),intent(in) :: mesh
+
+    if (associated(mesh%faces)) then
+      element_count=mesh%faces%unique_surface_element_count
+    else
+      element_count=0
+    end if
+
+  end function unique_surface_element_count_mesh  
+
   pure function face_count_scalar(field) result (face_count)
     ! Return the number of faces in a mesh.
     integer :: face_count

=== modified file 'femtools/Fields_Data_Types.F90'
--- femtools/Fields_Data_Types.F90	2011-10-06 15:16:29 +0000
+++ femtools/Fields_Data_Types.F90	2013-10-17 12:55:48 +0000
@@ -120,9 +120,27 @@
      integer, dimension(:), pointer :: coplanar_ids => null()
      !! a DG version of the surface mesh, useful for storing bc values
      type(mesh_type), pointer:: dg_surface_mesh => null()
-     !! A logical indicating if this mesh has internal boundaries
-     !! This means that element owners need to be written when writing out this mesh
-     logical :: has_internal_boundaries=.false.
+     !! A logical indicating if this mesh has a discontinuous internal boundary
+     !! This means that the pairs of internal facets are allowed to have two different
+     !! surface ids. When writing out this mesh both facets are written out and
+     !! element owners (an extra column indiciating which element is adjacent to each facet)
+     !! needs to be written out along with the surface mesh.
+     !! This is currently only used for periodic meshes which have a physical and aliased surface id
+     !! along the periodic boundary (which is an internal boundary of the periodic mesh).
+     !! Note that other meshes (with has_internal_boundaries==.false.) may still have internal facets
+     !! as part of the surface mesh, in this case the surface ids do have to agree and only one of each
+     !! pair of facets is written when writing out the mesh.
+     logical :: has_discontinuous_internal_boundaries=.false.
+     !! If internal facets are present in the surface mesh (and has_discontinuous_internal_boundaries==.false.) 
+     !! the surface facets are numbered such that 1:unique_surface_element_count visits each external facet
+     !! and each pair of internal facets only once (the order of this is typically determined by the 
+     !! read-in input mesh in which internal facets are only also only present once)
+     !! The second facets of each pair of internal facets are numbered 
+     !! unique_surface_element_count+1:surface_element_count (surface_element_count()==size(boundary_ids))
+     !! For meshes with no internal facets: unique_surface_element_count==surface_element_count()
+     !! For meshes with has_discontinuous_internal_boundaries no order is guaranteed and also
+     !! unique_surface_element_count==surface_element_count()
+     integer :: unique_surface_element_count
   end type mesh_faces
 
   type mesh_subdomain_mesh

=== modified file 'femtools/Surface_Labels.F90'
--- femtools/Surface_Labels.F90	2010-11-16 15:37:38 +0000
+++ femtools/Surface_Labels.F90	2013-10-17 12:55:48 +0000
@@ -282,74 +282,73 @@
   type(vector_field), intent(in):: positions
   integer, dimension(:), pointer:: coplanar_ids
     type(mesh_type), pointer:: surface_mesh
-    type(csr_sparsity), pointer :: eelist
     type(ilist) front
-    real, allocatable, dimension(:,:):: normals, normalgi
+    real, allocatable, dimension(:,:):: normalgi, normals
     real, allocatable, dimension(:):: detwei_f
-    integer, allocatable, dimension(:):: face_nodes
     real coplanar
-    integer, dimension(:), pointer:: neigh
+    integer, dimension(:), pointer:: neigh, nodes
     integer current_id, sngi, stotel
-    integer j, k, sele, pos, ele
+    integer j, k, sele, sele2, pos, ele
     
+    ewrite(1,*) "Inside get_coplanar_ids"
+
     if (.not. has_faces(mesh)) then
        call add_faces(mesh)
     end if
     
     stotel = surface_element_count(mesh)
+
+    if (.not. associated(mesh%faces%coplanar_ids)) then
+       allocate(mesh%faces%coplanar_ids(1:stotel))
+       coplanar_ids => mesh%faces%coplanar_ids
+    else
+       ! we assume they have been calculated already:
+       coplanar_ids => mesh%faces%coplanar_ids
+       return
+    end if
+
     if(stotel == 0) then
       sngi = 0
     else
       sngi = face_ngi(mesh, 1)
     end if
     surface_mesh => mesh%faces%surface_mesh
-    
-    if (.not. associated(mesh%faces%coplanar_ids)) then
-       allocate(mesh%faces%coplanar_ids(1:stotel))
-       coplanar_ids => mesh%faces%coplanar_ids
-    else
-       ! we assume they have been calculated already:
-       coplanar_ids => mesh%faces%coplanar_ids
-       return
-    end if
-    coplanar_ids => mesh%faces%coplanar_ids
-    allocate( normalgi(positions%dim,sngi), &
-       detwei_f(sngi), face_nodes(1:face_loc(mesh,1)) )
+
+    allocate(normalgi(positions%dim,sngi), detwei_f(sngi)) 
        
     ! Calculate element normals for all surface elements
     allocate(normals(positions%dim, stotel))
     do sele=1, stotel
        ele=face_ele(mesh, sele)
-       ! we don't want to use face_val on positions 
-       ! as it may not have a faces object (as opposed to the mesh)
-       face_nodes=face_global_nodes(mesh, sele)
        call transform_facet_to_physical( &
           positions, sele, detwei_f, normalgi)
        ! average over gauss points:
        normals(:,sele)=matmul(normalgi, detwei_f)/sum(detwei_f)
     end do
-    deallocate(normalgi, detwei_f, face_nodes)
+    deallocate(normalgi, detwei_f)
     
-    ! create element-element list for surface mesh
-    eelist => extract_eelist(surface_mesh)
+    ! create node-element list for surface mesh
+    ! (Note that we can't always construct an eelist, which would have been
+    ! more useful, because the surface mesh may split (for instance where an
+    ! external boundary meets an internal boundary) in which case multiple
+    ! surface elements can be connected via the same edge (facet of the surface element))
+    call add_nelist(surface_mesh)
     
     coplanar_ids = 0
     current_id = 1
     pos = 1
     do while (.true.)
-       ! Create a new starting point
+       ! Create a new starting point by finding a surface element without coplanar id
        do sele=pos, stotel
-          if(coplanar_ids(sele)==0) then
-             ! This is the first element in the new patch
-             pos = sele
-             coplanar_ids(pos) = current_id
-             exit
-          end if
+          if(coplanar_ids(sele)==0) exit
        end do
           
        ! Jump out of this while loop if we are finished
        if (sele>stotel) exit
           
+       ! This is the first element in the new patch
+       pos = sele
+       coplanar_ids(pos) = current_id
        ! Initialise the front
        call insert_ascending(front, pos)
           
@@ -358,19 +357,25 @@
           sele = pop(front)
           
           ! surrounding surface elements:
-          neigh => row_m_ptr(eelist, sele)
-          do j=1, size(neigh)
-             k = neigh(j)
-             if (k>0) then
-                if(coplanar_ids(k)==0) then
-                   coplanar = dot_product(normals(:,pos), normals(:,k))
+          ! note that we not only consider directly adjacent surface elements,
+          ! but also surface elements that only connect via one other node, as
+          ! long as they have a normal in the same direction however they should
+          ! be on the same plane
+          nodes => ele_nodes(surface_mesh, sele)
+          do j=1, size(nodes)
+             ! loop over surface elements connected to nodes(j)
+             neigh => node_neigh(surface_mesh, nodes(j))
+             do k=1, size(neigh)
+                sele2 = neigh(k)
+                if(coplanar_ids(sele2)==0) then
+                   coplanar = abs(dot_product(normals(:,pos), normals(:,sele2)))
                    if(coplanar>=COPLANAR_MAGIC_NUMBER) then
                    
-                      call insert_ascending(front, k)
-                      coplanar_ids(k) = current_id
+                      call insert_ascending(front, sele2)
+                      coplanar_ids(sele2) = current_id
                    end if
                 end if
-            end if
+             end do
           end do
        end do
          
@@ -379,6 +384,8 @@
     end do
     deallocate(normals)
     
+    ewrite(2,*) "Before merge_surface_ids, n/o local coplanes:", current_id-1
+
     call merge_surface_ids(mesh, coplanar_ids, max_id = current_id - 1)
 
   end subroutine get_coplanar_ids

=== modified file 'femtools/Write_GMSH.F90'
--- femtools/Write_GMSH.F90	2011-02-09 19:08:48 +0000
+++ femtools/Write_GMSH.F90	2013-10-17 12:55:48 +0000
@@ -258,12 +258,12 @@
     integer :: e, f, elemID
     character, parameter :: newLineChar=char(10)
 
-    logical :: internal_boundaries
+    logical :: needs_element_owners
 
     ! Gather some info about the mesh
     numElements = ele_count(mesh)
-    numFaces = surface_element_count(mesh)
-    internal_boundaries = has_internal_boundaries(mesh)
+    numFaces = unique_surface_element_count(mesh)
+    needs_element_owners = has_discontinuous_internal_boundaries(mesh)
 
     ! In the GMSH format, faces are also elements.
     numGMSHElems = numElements + numFaces
@@ -320,7 +320,7 @@
     ! Faces written out first
     
     ! Number of tags associated with elements
-    if(internal_boundaries) then
+    if(needs_element_owners) then
       ! write surface id and element owner
       numTags = 4
     else

=== modified file 'femtools/Write_Triangle.F90'
--- femtools/Write_Triangle.F90	2011-02-09 19:08:48 +0000
+++ femtools/Write_Triangle.F90	2013-10-17 12:55:48 +0000
@@ -233,18 +233,17 @@
     type(mesh_type), intent(in):: mesh
     
     integer unit, dim, nofaces, i
-    integer :: stotel, dg_total ! dg_total counts each internal face twice
+    integer :: dg_total ! dg_total counts each internal face twice
     integer :: ele, neigh, j, face
     integer, dimension(:), pointer :: neighbours
     
     unit=free_unit()
     
     dim=mesh_dim(mesh)
-    stotel=surface_element_count(mesh)
     dg_total=size(mesh%faces%face_list%sparsity%colm)
     
-    nofaces = (dg_total - stotel) / 2 ! internal faces, only once
-    nofaces = nofaces + stotel ! and the surface mesh
+    nofaces = (dg_total - surface_element_count(mesh)) / 2 ! internal faces, only once
+    nofaces = nofaces + unique_surface_element_count(mesh) ! and the surface mesh (again counting internal facets only once)
     
     select case(dim)
       case(3)
@@ -301,7 +300,7 @@
     unit=free_unit()
     
     dim=mesh_dim(mesh)
-    nofaces=surface_element_count(mesh)
+    nofaces=unique_surface_element_count(mesh)
     
     select case(dim)
       case(3)
@@ -315,7 +314,7 @@
         FLAbort("Invalid dimension")
     end select
 
-    if (has_internal_boundaries(mesh)) then
+    if (has_discontinuous_internal_boundaries(mesh)) then
       ! If the mesh is periodic, we want to write out the parent element of every face
       nolabels = 2
     else
@@ -325,7 +324,7 @@
     ! header line: nofaces, and number of boundary markers
     write(unit, *, err=42) nofaces, nolabels
     
-    if (.not. has_internal_boundaries(mesh)) then
+    if (.not. has_discontinuous_internal_boundaries(mesh)) then
       do i=1, nofaces
          write(unit, *, err=42) i, face_global_nodes(mesh, i), &
               surface_element_id(mesh, i)

=== modified file 'manual/configuring_fluidity.tex'
--- manual/configuring_fluidity.tex	2013-10-11 10:48:07 +0000
+++ manual/configuring_fluidity.tex	2013-10-17 12:55:48 +0000
@@ -1666,6 +1666,27 @@
 imposes a weak prescribed normal flow boundary condition on the surfaces specified.
 
 
+\subsection{Internal boundary conditions}\label{sec:internal_bcs}
+\index{boundary conditions!internal boundaries}
+Fluidity has some support (with limitations) for applying boundary conditions
+along internal boundaries. This is done in the normal way 
+by applying physical line/surface ids to the internal boundary, and associating
+boundary conditions with these in the options file. It should be
+noted that this functionality is not as well tested as other boundary
+conditions, so you should always first try it out on a simplified version of
+your problem to see if it behaves as expected. Strong Dirichlet boundary
+conditions should work correctly, other boundary condition types may not give
+the right answer. Since the gmsh .msh format does not give you the
+option of setting a different surface id on either side of the internal
+boundary, it is not possible to apply a different boundary condition type or
+value on either side.
+
+When using mesh adaptivity, it will maintain the internal boundary and its
+surface id. Mesh adaptivity in three dimensions with internal boundaries 
+is currently not supported. For parallel runs you have to use flredecomp (see
+section \ref{mesh!meshing tools!flredecomp}) to
+decompose the problem, as fldecomp will refuse any mesh with internal boundaries.
+
 \subsection{Special input date for boundary conditions}\label{sec:BCs:specialised}
 
 When running free surface simulations the surface elevation at the boundary is specified by applying a pressure Dirichlet condition. Since the free surface elevation is often measured data, there are some special possibilities to specify a pressure Dirichlet condition: 

=== modified file 'preprocessor/Populate_Sub_State.F90'
--- preprocessor/Populate_Sub_State.F90	2011-05-05 17:31:42 +0000
+++ preprocessor/Populate_Sub_State.F90	2013-10-17 12:55:48 +0000
@@ -67,7 +67,7 @@
 
   public populate_sub_state, set_full_domain_prescribed_fields, &
       &  sub_state_remap_to_full_mesh, update_subdomain_fields, &
-      &  use_sub_state
+      &  use_sub_state, populate_sub_state_module_check_options
 
 contains
 
@@ -205,6 +205,9 @@
       deallocate(prescribed_region_ids)
     end do
 
+    ! failing this may be caused by not preserving region ids in adaptivity - see options check below
+    assert(associated(external_mesh%region_ids))
+
     ! Derive subdomain_mesh_element list (an integer set):
     call allocate(subdomain_mesh_element_list)
     do ele = 1, element_count(external_mesh)
@@ -640,4 +643,19 @@
 
   end subroutine sub_state_remap_to_full_mesh
 
+  subroutine populate_sub_state_module_check_options
+
+    if (option_count('/material_phase/vector_field::Velocity/prognostic/prescribed_region')>0) then
+
+      if (have_option('/mesh_adaptivity/hr_adaptivity') .and. &
+        .not. have_option('/mesh_adaptivity/hr_adaptivity/preserve_mesh_regions')) then
+
+        ewrite(0,*) "When using prescribed regions with mesh adaptivity:"
+        FLExit("Need /mesh_adaptivity/hr_adaptivity/preserve_mesh_regions option")
+        
+      end if
+
+    end if
+  end subroutine populate_sub_state_module_check_options
+
 end module populate_sub_state_module

=== modified file 'tests/Stokes_Subduction_VK_Case_1b/Makefile'
--- tests/Stokes_Subduction_VK_Case_1b/Makefile	2011-10-21 10:05:25 +0000
+++ tests/Stokes_Subduction_VK_Case_1b/Makefile	2013-10-17 12:55:48 +0000
@@ -1,6 +1,5 @@
 input: clean
 	cp src/Subduction_Mesh.msh ./
-	../../bin/gmsh2triangle --2d -i Subduction_Mesh.msh
 
 clean:
 	rm -f *.ele *.edge *.node *.vtu *.stat *.msh *.detectors *.dat fluidity.* *_state matrixdump*

=== modified file 'tests/Stokes_Subduction_VK_Case_1b/Stokes-Subduction-VK-Case-1b.flml'
--- tests/Stokes_Subduction_VK_Case_1b/Stokes-Subduction-VK-Case-1b.flml	2013-01-06 19:51:01 +0000
+++ tests/Stokes_Subduction_VK_Case_1b/Stokes-Subduction-VK-Case-1b.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="Subduction_Mesh">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/flredecomp_2d_fieldweighted/Makefile'
--- tests/flredecomp_2d_fieldweighted/Makefile	2012-09-10 14:14:35 +0000
+++ tests/flredecomp_2d_fieldweighted/Makefile	2013-10-17 12:55:48 +0000
@@ -3,9 +3,8 @@
 default: input
 
 input:  clean
-	gmsh -2 -optimize src/Subduction_Mesh.geo
-	cp src/Subduction_Mesh.msh .
-	../../bin/gmsh2triangle --2d -i Subduction_Mesh.msh
+	gmsh -2 -optimize src/Subduction_Mesh.geo -o Subduction_Mesh.msh
+	../../bin/gmsh2triangle --2d Subduction_Mesh.msh
 
 clean:
 	rm -f  *.ele *.edge *.node *.vtu *.stat *.msh *.detectors fluidity.* Parallel*

=== modified file 'tests/flredecomp_2d_fieldweighted/flredecomp-2d-fieldweighted.flml'
--- tests/flredecomp_2d_fieldweighted/flredecomp-2d-fieldweighted.flml	2012-10-17 15:59:53 +0000
+++ tests/flredecomp_2d_fieldweighted/flredecomp-2d-fieldweighted.flml	2013-10-17 12:55:48 +0000
@@ -2,7 +2,7 @@
 <fluidity_options>
   <simulation_name>
     <string_value lines="1">flredecomp_field_weighted_partition_test</string_value>
-    <comment>Checks for field_weighted_partitions in flredecomp and gmsh2triangle -i.</comment>
+    <comment>Checks for field_weighted_partitions in flredecomp</comment>
   </simulation_name>
   <problem_type>
     <string_value lines="1">stokes</string_value>

=== added directory 'tests/internal_boundary_adaptivity'
=== added file 'tests/internal_boundary_adaptivity/Makefile'
--- tests/internal_boundary_adaptivity/Makefile	1970-01-01 00:00:00 +0000
+++ tests/internal_boundary_adaptivity/Makefile	2013-10-17 12:55:48 +0000
@@ -0,0 +1,6 @@
+input: clean
+	gmsh -2 src/square.geo -o square.msh
+
+clean:
+	rm -rf *.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele *.halo *.pvtu viscous_fs_zhong_B_* viscous_fs_zhong_A_* fluidity.* *.pyc
+	rm -rf *checkpoint* matrixdump*

=== added file 'tests/internal_boundary_adaptivity/internal_boundary.flml'
--- tests/internal_boundary_adaptivity/internal_boundary.flml	1970-01-01 00:00:00 +0000
+++ tests/internal_boundary_adaptivity/internal_boundary.flml	2013-10-17 12:55:48 +0000
@@ -0,0 +1,391 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">ib</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">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="DGMesh">
+      <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">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="CoordinateMesh"/>
+    <checkpointing>
+      <checkpoint_period_in_dumps>
+        <integer_value rank="0">1000</integer_value>
+      </checkpoint_period_in_dumps>
+      <checkpoint_at_end/>
+    </checkpointing>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.1</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">0.5</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="Fluid">
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_form/>
+            </stress_terms>
+          </continuous_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">1e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">200</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>
+        <prescribed_region name="Bottom">
+          <region_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </region_ids>
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </prescribed_region>
+        <boundary_conditions name="TopLeft">
+          <surface_ids>
+            <integer_value shape="1" rank="1">9</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">1.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">8</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <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="Scalar" rank="0">
+      <prognostic>
+        <mesh name="DGMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+            <diffusion_scheme>
+              <bassi_rebay/>
+            </diffusion_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>
+        </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">1e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">200</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Inflow">
+          <surface_ids>
+            <integer_value shape="1" rank="1">9</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Bottom">
+          <surface_ids>
+            <integer_value shape="1" rank="1">6</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="InternalBoundary">
+          <surface_ids>
+            <integer_value shape="1" rank="1">10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Diffusivity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.1</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="AnalyticalSolutionBottom" rank="0">
+      <prescribed>
+        <mesh name="DGMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value type="code" lines="20" language="python">def val(X,t):
+  return X[1]*2.0</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="BottomOnly" rank="0">
+      <prescribed>
+        <mesh name="DGMesh"/>
+        <value name="Top">
+          <region_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </region_ids>
+          <constant>
+            <real_value rank="0">0.0</real_value>
+          </constant>
+        </value>
+        <value name="Bottom">
+          <region_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </region_ids>
+          <constant>
+            <real_value rank="0">1.0</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="Error" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value type="code" lines="20" language="python">asb = state.scalar_fields['AnalyticalSolutionBottom']
+bo = state.scalar_fields['BottomOnly']
+s = state.scalar_fields['Scalar']
+field.val[:] = bo.val[:]*(asb.val[:]-s.val[:])</string_value>
+          <depends>
+            <string_value lines="1">AnalyticalSolutionBottom,BottomOnly</string_value>
+          </depends>
+        </algorithm>
+        <mesh name="DGMesh"/>
+        <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>
+  <mesh_adaptivity>
+    <mesh_movement>
+      <imposed_grid_velocity/>
+      <vector_field name="GridVelocity" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import sin,cos,pi
+  v = 0.1*sin(2*pi*X[0])*sin(2*pi*t)*cos(2*pi*(X[1]-0.5))
+  return 0.0, v</string_value>
+            </python>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </mesh_movement>
+    <hr_adaptivity>
+      <period_in_timesteps>
+        <integer_value rank="0">3</integer_value>
+      </period_in_timesteps>
+      <maximum_number_of_nodes>
+        <integer_value rank="0">100000</integer_value>
+      </maximum_number_of_nodes>
+      <enable_gradation/>
+      <tensor_field name="MinimumEdgeLengths">
+        <anisotropic_symmetric>
+          <constant>
+            <real_value symmetric="true" dim2="dim" shape="2 2" dim1="dim" rank="2">0.1 0.0 0.0 0.1</real_value>
+          </constant>
+        </anisotropic_symmetric>
+      </tensor_field>
+      <tensor_field name="MaximumEdgeLengths">
+        <anisotropic_symmetric>
+          <constant>
+            <real_value symmetric="true" dim2="dim" shape="2 2" dim1="dim" rank="2">0.1 0.0 0.0 0.1</real_value>
+          </constant>
+        </anisotropic_symmetric>
+      </tensor_field>
+      <preserve_mesh_regions/>
+    </hr_adaptivity>
+  </mesh_adaptivity>
+</fluidity_options>

=== added file 'tests/internal_boundary_adaptivity/internal_boundary_adaptivity.xml'
--- tests/internal_boundary_adaptivity/internal_boundary_adaptivity.xml	1970-01-01 00:00:00 +0000
+++ tests/internal_boundary_adaptivity/internal_boundary_adaptivity.xml	2013-10-17 12:55:48 +0000
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>interal_boundary_adaptivity<comment>Tests: internal boundaries+adaptivity+checkpointing</comment></name>
+  <owner userid="skramer"/>
+  <problem_definition length="short" nprocs="1">
+    <command_line>fluidity -v2 -l internal_boundary.flml &amp;&amp;
+spud-set ib_1_checkpoint.flml /timestepping/finish_time 0.999999 &amp;&amp;
+fluidity -v2 -l ib_1_checkpoint.flml</command_line>
+  </problem_definition>
+  <variables>
+    <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="stats" language="python">from fluidity_tools import stat_parser
+s=stat_parser("ib.stat")
+stats={}
+stats["nodes_begin"]=s["CoordinateMesh"]["nodes"][0]
+stats["nodes_end"]=s["CoordinateMesh"]["nodes"][-1]
+stats["asbmin"]=s["Fluid"]["AnalyticalSolutionBottom"]["min"][-1]
+stats["asbmax"]=s["Fluid"]["AnalyticalSolutionBottom"]["max"][-1]
+stats["vell2_begin"]=s["Fluid"]["Velocity%magnitude"]["l2norm"][0]
+stats["vell2_end"]=s["Fluid"]["Velocity%magnitude"]["l2norm"][-1]
+stats["scalarl2_begin"]=s["Fluid"]["Scalar"]["l2norm"][0]
+stats["scalarl2_end"]=s["Fluid"]["Scalar"]["l2norm"][-1]</variable>
+    <variable name="stats_after_checkpoint" language="python">from fluidity_tools import stat_parser
+s=stat_parser("ib_checkpoint.stat")
+stats_after_checkpoint={}
+stats_after_checkpoint["nodes_begin"]=s["CoordinateMesh"]["nodes"][0]
+stats_after_checkpoint["nodes_end"]=s["CoordinateMesh"]["nodes"][-1]
+stats_after_checkpoint["asbmin"]=s["Fluid"]["AnalyticalSolutionBottom"]["min"][-1]
+stats_after_checkpoint["asbmax"]=s["Fluid"]["AnalyticalSolutionBottom"]["max"][-1]
+stats_after_checkpoint["vell2_begin"]=s["Fluid"]["Velocity%magnitude"]["l2norm"][0]
+stats_after_checkpoint["vell2_end"]=s["Fluid"]["Velocity%magnitude"]["l2norm"][-1]
+stats_after_checkpoint["scalarl2_begin"]=s["Fluid"]["Scalar"]["l2norm"][0]
+stats_after_checkpoint["scalarl2_end"]=s["Fluid"]["Scalar"]["l2norm"][-1]
+stats_after_checkpoint["error"]=s["Fluid"]["Error"]["l2norm"][-1]</variable>
+  </variables>
+  <pass_tests>
+    <test name="Solvers converged" language="python">assert(solvers_converged)</test>
+    <test name="NodesDecreaseInAdapt" language="python">assert stats['nodes_begin']&gt;stats['nodes_end']</test>
+    <test name="ScalarL2Increased" language="python">assert stats['scalarl2_begin']&lt;stats['scalarl2_end']</test>
+    <test name="VelocityL2Increased" language="python">assert stats['vell2_begin']&lt;stats['vell2_end']</test>
+    <test name="NodesDecreaseInAdapt" language="python">assert stats['nodes_begin']&gt;stats['nodes_end']</test>
+    <test name="MeshMovementMax" language="python">assert abs(stats['asbmax']-2.06)&lt;1e-3</test>
+    <test name="SameNumberOfNodes" language="python">assert stats["nodes_end"]==stats_after_checkpoint["nodes_begin"]</test>
+    <test name="NearlySameScalarL2" language="python">assert abs(stats["scalarl2_end"]-stats_after_checkpoint["scalarl2_begin"])&lt;0.03</test>
+    <test name="CloseVelocityL2" language="python">assert abs(stats["vell2_end"]-stats_after_checkpoint["vell2_begin"])&lt;0.1</test>
+    <test name="ScalarL2IncreasedAfter" language="python">assert stats_after_checkpoint['scalarl2_begin']&lt;stats_after_checkpoint['scalarl2_end']</test>
+    <test name="NodesConstantAfter" language="python">assert abs(stats_after_checkpoint['nodes_begin']-stats_after_checkpoint['nodes_end'])&lt;10</test>
+    <test name="MeshMovementMin" language="python">assert abs(stats['asbmin']+0.06)&lt;1e-3</test>
+    <test name="MeshMovedBackTop" language="python">assert abs(stats_after_checkpoint['asbmax']-2.00)&lt;3e-3</test>
+    <test name="MeshMovedBackBottom" language="python">assert abs(stats_after_checkpoint['asbmin'])&lt;3e-3</test>
+    <test name="SmallError" language="python">assert abs(stats_after_checkpoint['error'])&lt;0.02</test>
+  </pass_tests>
+</testproblem>

=== added directory 'tests/internal_boundary_adaptivity/src'
=== added file 'tests/internal_boundary_adaptivity/src/square.geo'
--- tests/internal_boundary_adaptivity/src/square.geo	1970-01-01 00:00:00 +0000
+++ tests/internal_boundary_adaptivity/src/square.geo	2013-10-17 12:55:48 +0000
@@ -0,0 +1,33 @@
+l=20;
+p=1.0;
+Point(1) = {0, 0, 0};
+Extrude {1, 0, 0} {
+  Point{1}; 
+}
+Extrude {0, 0.5, 0} {
+  Line{1};
+}
+Transfinite Line{1} = l+1;
+Transfinite Line{2} = l+1;
+Transfinite Line{-3} = (l+2)/2 Using Progression p;
+Transfinite Line{-4} = (l+2)/2 Using Progression p;
+Transfinite Surface{5} = {1,2,3,4} Right;
+Extrude {0, 0.5, 0} {
+  Line{2};
+}
+Transfinite Line{6} = l+1;
+Transfinite Line{7} = (l+2)/2 Using Progression p;
+Transfinite Line{8} = (l+2)/2 Using Progression p;
+Transfinite Surface{9} = {3,4,5,6} Right;
+// Bottom
+Physical Line(6) = {1};
+Physical Line(7) = {4, 8};
+// Top
+Physical Line(8) = {6};
+// Left
+Physical Line(9) = {7}; // Top-Left
+Physical Line(11) = {3}; // Bottom-Left
+// Internal
+Physical Line(10) = {2};
+Physical Surface(1) = {9}; // Top
+Physical Surface(2) = {5}; // Bottom

=== modified file 'tests/viscous_fs_zhong_spatial/Makefile'
--- tests/viscous_fs_zhong_spatial/Makefile	2013-09-17 15:06:47 +0000
+++ tests/viscous_fs_zhong_spatial/Makefile	2013-10-17 12:55:48 +0000
@@ -1,8 +1,13 @@
 input: clean
+<<<<<<< TREE
 	gmsh -2 src/squareA.geo
 	../../bin/gmsh2triangle --2d -i src/squareA.msh
 	gmsh -2 src/squareB.geo
 	../../bin/gmsh2triangle --2d -i src/squareB.msh
+=======
+	gmsh -2 src/squareA.geo -o squareA.msh
+	gmsh -2 src/squareB.geo -o squareB.msh
+>>>>>>> MERGE-SOURCE
 
 clean:
-	rm -rf src/*.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele *.halo *.pvtu viscous_fs_zhong_B_* viscous_fs_zhong_A_* fluidity.* *.pyc
+	rm -rf *.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele *.halo *.pvtu viscous_fs_zhong_B_* viscous_fs_zhong_A_* fluidity.* *.pyc

=== modified file 'tests/viscous_fs_zhong_spatial/viscous_fs_zhong_A.flml'
--- tests/viscous_fs_zhong_spatial/viscous_fs_zhong_A.flml	2012-03-16 15:27:11 +0000
+++ tests/viscous_fs_zhong_spatial/viscous_fs_zhong_A.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareA">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial/viscous_fs_zhong_B.flml'
--- tests/viscous_fs_zhong_spatial/viscous_fs_zhong_B.flml	2012-03-16 15:27:11 +0000
+++ tests/viscous_fs_zhong_spatial/viscous_fs_zhong_B.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareB">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial_explicit/Makefile'
--- tests/viscous_fs_zhong_spatial_explicit/Makefile	2013-09-17 15:06:47 +0000
+++ tests/viscous_fs_zhong_spatial_explicit/Makefile	2013-10-17 12:55:48 +0000
@@ -1,8 +1,13 @@
 input: clean
+<<<<<<< TREE
 	gmsh -2 src/squareA.geo
 	../../bin/gmsh2triangle --2d -i src/squareA.msh
 	gmsh -2 src/squareB.geo
 	../../bin/gmsh2triangle --2d -i src/squareB.msh
+=======
+	gmsh -2 src/squareA.geo -o squareA.msh
+	gmsh -2 src/squareB.geo -o squareB.msh
+>>>>>>> MERGE-SOURCE
 
 clean:
-	rm -rf src/*.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele
+	rm -rf *.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele

=== modified file 'tests/viscous_fs_zhong_spatial_explicit/viscous_fs_zhong_A.flml'
--- tests/viscous_fs_zhong_spatial_explicit/viscous_fs_zhong_A.flml	2012-03-16 15:53:58 +0000
+++ tests/viscous_fs_zhong_spatial_explicit/viscous_fs_zhong_A.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareA">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial_explicit/viscous_fs_zhong_B.flml'
--- tests/viscous_fs_zhong_spatial_explicit/viscous_fs_zhong_B.flml	2012-03-16 15:53:58 +0000
+++ tests/viscous_fs_zhong_spatial_explicit/viscous_fs_zhong_B.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareB">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial_explicit_varrho/Makefile'
--- tests/viscous_fs_zhong_spatial_explicit_varrho/Makefile	2013-09-17 15:06:47 +0000
+++ tests/viscous_fs_zhong_spatial_explicit_varrho/Makefile	2013-10-17 12:55:48 +0000
@@ -1,8 +1,13 @@
 input: clean
+<<<<<<< TREE
 	gmsh -2 src/squareA.geo
 	../../bin/gmsh2triangle --2d -i src/squareA.msh
 	gmsh -2 src/squareB.geo
 	../../bin/gmsh2triangle --2d -i src/squareB.msh
+=======
+	gmsh -2 src/squareA.geo -o squareA.msh
+	gmsh -2 src/squareB.geo -o squareB.msh
+>>>>>>> MERGE-SOURCE
 
 clean:
-	rm -rf src/*.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele
+	rm -rf *.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele

=== modified file 'tests/viscous_fs_zhong_spatial_explicit_varrho/viscous_fs_zhong_A.flml'
--- tests/viscous_fs_zhong_spatial_explicit_varrho/viscous_fs_zhong_A.flml	2012-03-16 05:00:48 +0000
+++ tests/viscous_fs_zhong_spatial_explicit_varrho/viscous_fs_zhong_A.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareA">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial_explicit_varrho/viscous_fs_zhong_B.flml'
--- tests/viscous_fs_zhong_spatial_explicit_varrho/viscous_fs_zhong_B.flml	2012-03-16 05:00:48 +0000
+++ tests/viscous_fs_zhong_spatial_explicit_varrho/viscous_fs_zhong_B.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareB">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial_varrho/Makefile'
--- tests/viscous_fs_zhong_spatial_varrho/Makefile	2013-09-17 15:02:30 +0000
+++ tests/viscous_fs_zhong_spatial_varrho/Makefile	2013-10-17 12:55:48 +0000
@@ -1,8 +1,13 @@
 input: clean
+<<<<<<< TREE
 	gmsh -2 src/squareA.geo
 	../../bin/gmsh2triangle --2d -i src/squareA.msh
 	gmsh -2 src/squareB.geo
 	../../bin/gmsh2triangle --2d -i src/squareB.msh
+=======
+	gmsh -2 src/squareA.geo -o squareA.msh
+	gmsh -2 src/squareB.geo -o squareB.msh
+>>>>>>> MERGE-SOURCE
 
 clean:
-	rm -rf src/*.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele
+	rm -rf *.msh *.stat *.vtu *.log-0 *.err-0 *checkpoint *.convergence *.log *.node *.edge *.ele

=== modified file 'tests/viscous_fs_zhong_spatial_varrho/viscous_fs_zhong_A.flml'
--- tests/viscous_fs_zhong_spatial_varrho/viscous_fs_zhong_A.flml	2012-03-16 15:53:58 +0000
+++ tests/viscous_fs_zhong_spatial_varrho/viscous_fs_zhong_A.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareA">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>

=== modified file 'tests/viscous_fs_zhong_spatial_varrho/viscous_fs_zhong_B.flml'
--- tests/viscous_fs_zhong_spatial_varrho/viscous_fs_zhong_B.flml	2012-03-16 15:53:58 +0000
+++ tests/viscous_fs_zhong_spatial_varrho/viscous_fs_zhong_B.flml	2013-10-17 12:55:48 +0000
@@ -12,7 +12,7 @@
     </dimension>
     <mesh name="CoordinateMesh">
       <from_file file_name="squareB">
-        <format name="triangle"/>
+        <format name="gmsh"/>
         <stat>
           <include_in_stat/>
         </stat>