Incorrect #ifndef in find_adjacent_element

[Mark Filipiak]

In find_adjacent_element() (femtools/Adjacency_Lists.F90, line 911) there is a #ifndef which (probably?) should be an #ifdef:

      adj_ele = 0
      ele_loop: do i = 1, size(elements1)
         candidate_ele = elements1(i)
         if(candidate_ele == ele) cycle ele_loop ! Ignore the query element
         ! See if this element borders all other nodes
         do j = 2, size(nodes)
           ! If candidate_ele is not in all row_idx, nodes are not boundary
           ! nodes for candidate_ele, and it isn't the adjacent element.
           if(.not. any(row_idx(j - 1)%ptr == candidate_ele)) cycle ele_loop
         end do

#ifndef DDEBUG <<------------------------------------------------------------------------------------------------------------------------------------------------
         if(adj_ele > 0) then
           ! We've found more than one adjacent element. We're in trouble.
           adj_ele = -1
           return
         end if
#endif
         adj_ele = candidate_ele
#ifndef DDEBUG
         ! We've found the adjacent element. We're done.
         return
#endif
      end do ele_loop

otherwise there is no test for more than one adjacent element with --enable-debugging, nor without --enable-debugging, since the second (correct) #ifndef exits the loop when the first candidate_ele is found.

[Mark Filipiak (mjf-q)]

And when the #ifndef is changed to #ifdef the following tests fail with --enable-debugging:

periodic_adaptivity.xml: PFFP
periodic_adaptivity_3d.xml: PFFP
periodic_2plus1_adaptivity.xml: PFFP
Stokes_Subduction_VK_Case_1b.xml: F
periodic_adaptivity_tracer.xml: PFFP

[Mark Filipiak (mjf-q)]

The error in Stokes_Subduction_VK_Case_1b is for the surfacemesh (i.e. line mesh since this is a 2d problem) and in the .edge file there are 4 edges that meet at node 3 , one edge for surface 11, one for 12, and two for 32 (and these seem to be a double edge!)
Surfaces meeting should be OK but will show more than 2 elements with one face (as for some periodic mesh cases) which will cause Fluidity to fail.

[Mark Filipiak (mjf-q)]

That double edge is correct, it is an internal surface.

[Stephan Kramer (s-kramer)]

Ah good catch. That probably also explains the periodic cases indeed, they have an internal boundary as well. We should probably look into why we call add_faces() on these surface_meshes (which I presume is where this shows up in). I can imagine it needing things like ele_neigh() on a surface mesh for some algorithms, but in these cases here it may be computing it incorrectly...

[Mark Filipiak (mjf-q)]

In the hilbert branch, I have already replaced the add_faces() call in renumber_positions_elements (Fields_Manipulation.F90) with a manual copy of the face list because the same problem occurs for periodic meshes with only 2 elements in the periodic direction.

[Mark Filipiak (mjf-q)]

Also tests

saltfinger_adapt_periodic.xml: PF
flredecomp_2d_fieldweighted.xml: F

[Stephan Kramer (s-kramer)]

I don't think we should be supporting periodic meshes with only 2 elements across. In 2d, its surface mesh would only have 2 nodes and two elements that are in fact exactly the same. This is quite unlikely to work, and would probably break most sanity checks.

[Mark Filipiak (mjf-q)]

There is one periodic test case with 2 elements across:

ekman.xml

[Mark Filipiak (mjf-q)]

And the tests recently added to the trunk:

viscous_fs_zhong_spatial.xml: F (this fails in flredecomp)
viscous_fs_zhong_spatial_varrho.xml: F
viscous_fs_zhong_spatial_explicit_varrho.xml: F
viscous_fs_zhong_spatial_explicit.xml: F

[Mark Filipiak (mjf-q)]

I meant that the viscous_fs_zhong_spatial* tests have the find_adjacent_element() problem, not that they have 2 elements in the periodic direction.

[Cian Wilson (cwilson)]

The zhong cases also have an internal boundary with a boundary id.