Merge lp:~wence/fluidity/fix-mom-dg-assembly into lp:fluidity

[Lawrence Mitchell]

I think this is the necessary fix for Christian's issue. Have not checked tests, so a buildbot branch for that would be good.

[Lawrence Mitchell (wence)]

Buildbot pass: http://buildbot-ocean.ese.ic.ac.uk:8080/builders/fix-mom-dg-assembly

[Stephan Kramer (s-kramer)]

That looks good. Your explanation in the comment makes sense to me.

review: Approve

=== modified file 'assemble/Momentum_DG.F90'
--- assemble/Momentum_DG.F90	2012-04-17 11:31:55 +0000
+++ assemble/Momentum_DG.F90	2012-08-13 12:49:20 +0000
@@ -905,7 +905,11 @@
     
     ! In parallel, we construct terms on elements we own and those in
     ! the L1 element halo.
-    assemble_element = .not.dg.or.element_neighbour_owned(U, ele)
+    ! Note that element_neighbour_owned(U, ele) may return .false. if
+    ! ele is owned.  For example, if ele is the only owned element on
+    ! this process.  Hence we have to check for element ownership
+    ! directly as well.
+    assemble_element = .not.dg.or.element_neighbour_owned(U, ele).or.element_owned(U, ele)
 
     primal = .not.dg
     if(viscosity_scheme == CDG) primal = .true.

=== modified file 'femtools/Parallel_fields.F90'
--- femtools/Parallel_fields.F90	2011-07-28 12:16:05 +0000
+++ femtools/Parallel_fields.F90	2012-08-13 12:49:20 +0000
@@ -284,8 +284,11 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in MESH that
     !!< is owned by this process otherwise .false.
 
-    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
-    !! halo.
+    !! Note, you cannot use this function to compute if
+    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
+    !! element on a process, then none of the neighbours will be owned.
+    !! You can use this function to compute whether ELEMENT_NUMBER is
+    !! in the L1 element halo.
     type(mesh_type), intent(in) :: mesh
     integer, intent(in) :: element_number
     logical :: owned
@@ -311,8 +314,11 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in FIELD that
     !!< is owned by this process otherwise .false.
 
-    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
-    !! halo.
+    !! Note, you cannot use this function to compute if
+    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
+    !! element on a process, then none of the neighbours will be owned.
+    !! You can use this function to compute whether ELEMENT_NUMBER is
+    !! in the L1 element halo.
     type(scalar_field), intent(in) :: field
     integer, intent(in) :: element_number
     logical :: owned
@@ -324,8 +330,11 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in FIELD that
     !!< is owned by this process otherwise .false.
 
-    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
-    !! halo.
+    !! Note, you cannot use this function to compute if
+    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
+    !! element on a process, then none of the neighbours will be owned.
+    !! You can use this function to compute whether ELEMENT_NUMBER is
+    !! in the L1 element halo.
     type(vector_field), intent(in) :: field
     integer, intent(in) :: element_number
     logical :: owned
@@ -337,8 +346,11 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in FIELD that
     !!< is owned by this process otherwise .false.
 
-    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
-    !! halo.
+    !! Note, you cannot use this function to compute if
+    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
+    !! element on a process, then none of the neighbours will be owned.
+    !! You can use this function to compute whether ELEMENT_NUMBER is
+    !! in the L1 element halo.
     type(tensor_field), intent(in) :: field
     integer, intent(in) :: element_number
     logical :: owned