quadtree.cc

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Non-inline and non-templated functions for QuadTree and QuadTreeForest
// classes
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include <set>
 
// oomph-lib headers
#include "quadtree.h"
#include "refineable_quad_element.h"
 
 
namespace oomph
{
  //========================================================================
  /// Bool indicating that static member data has been setup
  //========================================================================
  bool QuadTree::Static_data_has_been_setup = false;
 
  //========================================================================
  /// Colours for neighbours in various directions (static data)
  //========================================================================
  Vector<std::string> QuadTree::Colour;
 
  //=======================================================================
  ///  S_base(i,direction):  Initial value for coordinate s[i] on
  /// the edge indicated by direction (S/E/N/W) (static data)
  //======================================================================
  DenseMatrix<double> QuadTree::S_base;
 
  //======================================================================
  /// S_step(i,direction) Increments for coordinate s[i] when
  /// progressing along the edge indicated by direction (S/E/N/W);
  /// Left/lower vertex: S_base; Right/upper vertex: S_base + S_step
  /// (static data)
  //=====================================================================
  DenseMatrix<double> QuadTree::S_step;
 
  //=====================================================================
  /// Translate (enumerated) directions into strings (static data)
  //=====================================================================
  Vector<std::string> QuadTree::Direct_string;
 
  //====================================================================
  /// Get opposite edge, e.g. Reflect_edge[N]=S (static data)
  //====================================================================
  Vector<int> QuadTree::Reflect_edge;
 
  //===================================================================
  /// Array of direction/quadrant adjacency scheme:
  /// Is_adjacent(i_vertex_or_edge,j_quadrant): Is edge/vertex
  /// adjacent to quadrant? (static data)
  //===================================================================
  DenseMatrix<bool> QuadTree::Is_adjacent;
 
  //====================================================================
  /// Reflection scheme: Reflect(direction,quadrant): Get mirror
  /// of quadrant in specified direction. E.g. Reflect(S,NE)=SE
  /// (static data)
  //====================================================================
  DenseMatrix<int> QuadTree::Reflect;
 
  //=====================================================================
  /// Rotate coordinates: If north becomes NorthIs then direction
  /// becomes Rotate(NorthIs,direction). E.g.  Rotate(E,NW)=NE;
  /// (static data)
  //=====================================================================
  DenseMatrix<int> QuadTree::Rotate;
 
  //=====================================================================
  /// Angle betwen rotated coordinates: If old_direction becomes
  /// new_direction then the angle between the axes (in anti-clockwise
  /// direction is Rotate_angle(old_direction,new_direction); E.g.
  /// Rotate_angle(E,N)=90; (static data)
  //=====================================================================
  DenseMatrix<int> QuadTree::Rotate_angle;
 
  //==========================================================================
  /// S_direct(direction,son_quadrant): The lower left corner
  /// of son_quadrant has an offset of h/2 S_direct(direction,son_quadrant)
  /// in the specified direction. E.g. S_direct(S,NE)=1 and  S_direct(S,NW)=0
  /// (static data)
  //==========================================================================
  DenseMatrix<int> QuadTree::S_direct;
 
 
  //====================================================================
  /// Setup the static data stored in the QuadTree -- this needs to be
  /// called before QuadTrees can be used. Automatically called
  /// by RefineableQuadMesh constructor.
  //====================================================================
  void QuadTree::setup_static_data()
  {
    using namespace QuadTreeNames;
 
 
#ifdef PARANOID
    if (Tree::OMEGA != QuadTree::OMEGA)
    {
      std::ostringstream error_stream;
      error_stream << "Inconsistent enumeration!  \n    Tree::OMEGA="
                   << Tree::OMEGA << "\nQuadTree::OMEGA=" << QuadTree::OMEGA
                   << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
 
    // Set flag to indicate that static data has been setup
    Static_data_has_been_setup = true;
 
    // Tecplot Colours for neighbours in various directions
    Colour.resize(27);
    Colour[SW] = "YELLOW";
    Colour[SE] = "YELLOW";
    Colour[NW] = "YELLOW";
    Colour[NE] = "YELLOW";
    Colour[E] = "CYAN";
    Colour[W] = "RED";
    Colour[N] = "GREEN";
    Colour[S] = "BLUE";
    Colour[OMEGA] = "YELLOW";
 
    // S_base(i,direction):  Initial value for coordinate s[i] on the
    // edge indicated by direction (S/E/N/W)
    //
    // S_step(i,direction) Increments for coordinate s[i] when progressing
    // along the edge indicated by direction (S/E/N/W)
    S_base.resize(2, 27);
    S_step.resize(2, 27);
 
    S_base(0, N) = -1.0;
    S_base(1, N) = 1.0;
    S_step(0, N) = 2.0;
    S_step(1, N) = 0.0;
 
    S_base(0, S) = -1.0;
    S_base(1, S) = -1.0;
    S_step(0, S) = 2.0;
    S_step(1, S) = 0.0;
 
    S_base(0, W) = -1.0;
    S_base(1, W) = -1.0;
    S_step(0, W) = 0.0;
    S_step(1, W) = 2.0;
 
    S_base(0, E) = 1.0;
    S_base(1, E) = -1.0;
    S_step(0, E) = 0.0;
    S_step(1, E) = 2.0;
 
    // Translate (enumerated) directions into strings
    Direct_string.resize(27);
 
    Direct_string[SW] = "SW";
    Direct_string[NW] = "NW";
    Direct_string[SE] = "SE";
    Direct_string[NE] = "NE";
 
    Direct_string[S] = "S";
    Direct_string[N] = "N";
    Direct_string[E] = "E";
    Direct_string[W] = "W";
    Direct_string[OMEGA] = "OMEGA";
 
    // Build direction/quadrant adjacency scheme
    // Is_adjacent(i_vertex_or_edge,j_quadrant):
    // Is edge/vertex adjacent to quadrant?
    Is_adjacent.resize(27, 27);
 
    Is_adjacent(N, NW) = true;
    Is_adjacent(E, NW) = false;
    Is_adjacent(S, NW) = false;
    Is_adjacent(W, NW) = true;
    Is_adjacent(NW, NW) = true;
    Is_adjacent(NE, NW) = false;
    Is_adjacent(SW, NW) = false;
    Is_adjacent(SE, NW) = false;
 
    Is_adjacent(N, NE) = true;
    Is_adjacent(E, NE) = true;
    Is_adjacent(S, NE) = false;
    Is_adjacent(W, NE) = false;
    Is_adjacent(NW, NE) = false;
    Is_adjacent(NE, NE) = true;
    Is_adjacent(SW, NE) = false;
    Is_adjacent(SE, NE) = false;
 
    Is_adjacent(N, SW) = false;
    Is_adjacent(E, SW) = false;
    Is_adjacent(S, SW) = true;
    Is_adjacent(W, SW) = true;
    Is_adjacent(NW, SW) = false;
    Is_adjacent(NE, SW) = false;
    Is_adjacent(SW, SW) = true;
    Is_adjacent(SE, SW) = false;
 
 
    Is_adjacent(N, SE) = false;
    Is_adjacent(E, SE) = true;
    Is_adjacent(S, SE) = true;
    Is_adjacent(W, SE) = false;
    Is_adjacent(NW, SE) = false;
    Is_adjacent(NE, SE) = false;
    Is_adjacent(SW, SE) = false;
    Is_adjacent(SE, SE) = true;
 
    // Rotation scheme: If north becomes NorthIs then
    // direction becomes Rotate(NorthIs,direction)
    // Initialise to OMEGA
    Rotate.resize(27, 27);
 
    Rotate(N, N) = N;
    Rotate(N, E) = E;
    Rotate(N, S) = S;
    Rotate(N, W) = W;
    Rotate(N, NW) = NW;
    Rotate(N, NE) = NE;
    Rotate(N, SE) = SE;
    Rotate(N, SW) = SW;
 
    Rotate(W, N) = W;
    Rotate(W, E) = N;
    Rotate(W, S) = E;
    Rotate(W, W) = S;
    Rotate(W, NW) = SW;
    Rotate(W, NE) = NW;
    Rotate(W, SE) = NE;
    Rotate(W, SW) = SE;
 
    Rotate(S, N) = S;
    Rotate(S, E) = W;
    Rotate(S, S) = N;
    Rotate(S, W) = E;
    Rotate(S, NW) = SE;
    Rotate(S, NE) = SW;
    Rotate(S, SE) = NW;
    Rotate(S, SW) = NE;
 
    Rotate(E, N) = E;
    Rotate(E, E) = S;
    Rotate(E, S) = W;
    Rotate(E, W) = N;
    Rotate(E, NW) = NE;
    Rotate(E, NE) = SE;
    Rotate(E, SE) = SW;
    Rotate(E, SW) = NW;
 
    // Angle betwen rotated coordinates:
    // old_direction becomes new_direction then the angle between
    // the axes is Rotate_angle(old_direction,new_direction)
    Rotate_angle.resize(27, 27);
 
    Rotate_angle(N, N) = 0;
    Rotate_angle(N, W) = 90;
    Rotate_angle(N, S) = 180;
    Rotate_angle(N, E) = 270;
 
    Rotate_angle(S, S) = 0;
    Rotate_angle(S, E) = 90;
    Rotate_angle(S, N) = 180;
    Rotate_angle(S, W) = 270;
 
    Rotate_angle(W, W) = 0;
    Rotate_angle(W, S) = 90;
    Rotate_angle(W, E) = 180;
    Rotate_angle(W, N) = 270;
 
    Rotate_angle(E, E) = 0;
    Rotate_angle(E, N) = 90;
    Rotate_angle(E, W) = 180;
    Rotate_angle(E, S) = 270;
 
    // Reflection scheme:
    // Reflect(direction,quadrant): Get mirror of quadrant in direction
    Reflect.resize(27, 27);
 
    Reflect(N, NW) = SW;
    Reflect(E, NW) = NE;
    Reflect(S, NW) = SW;
    Reflect(W, NW) = NE;
    Reflect(NW, NW) = SE;
    Reflect(NE, NW) = SE;
    Reflect(SW, NW) = SE;
    Reflect(SE, NW) = SE;
 
    Reflect(N, NE) = SE;
    Reflect(E, NE) = NW;
    Reflect(S, NE) = SE;
    Reflect(W, NE) = NW;
    Reflect(NW, NE) = SW;
    Reflect(NE, NE) = SW;
    Reflect(SW, NE) = SW;
    Reflect(SE, NE) = SW;
 
    Reflect(N, SW) = NW;
    Reflect(E, SW) = SE;
    Reflect(S, SW) = NW;
    Reflect(W, SW) = SE;
    Reflect(NW, SW) = NE;
    Reflect(NE, SW) = NE;
    Reflect(SW, SW) = NE;
    Reflect(SE, SW) = NE;
 
    Reflect(N, SE) = NE;
    Reflect(E, SE) = SW;
    Reflect(S, SE) = NE;
    Reflect(W, SE) = SW;
    Reflect(NW, SE) = NW;
    Reflect(NE, SE) = NW;
    Reflect(SW, SE) = NW;
    Reflect(SE, SE) = NW;
 
    // S_direct(direction,son_quadrant): The lower left corner of
    // my son_quadrant has an offset of h/2 S_direct(direction,son_quadrant)
    // in the direction. [Look in the positive direction along the
    // edge 'direction' and measure the offset of the lower left corner
    // of the son quadrant in this direction]
    S_direct.resize(27, 27);
 
    S_direct(W, NW) = 1;
    S_direct(E, NW) = 1;
    S_direct(N, NW) = 0;
    S_direct(S, NW) = 0;
 
    S_direct(W, NE) = 1;
    S_direct(E, NE) = 1;
    S_direct(N, NE) = 1;
    S_direct(S, NE) = 1;
 
    S_direct(W, SW) = 0;
    S_direct(E, SW) = 0;
    S_direct(N, SW) = 0;
    S_direct(S, SW) = 0;
 
    S_direct(W, SE) = 0;
    S_direct(E, SE) = 0;
    S_direct(N, SE) = 1;
    S_direct(S, SE) = 1;
 
    // Get opposite edge, e.g. Reflect_edge(N)=S
    Reflect_edge.resize(27);
 
    Reflect_edge[N] = S;
    Reflect_edge[S] = N;
    Reflect_edge[E] = W;
    Reflect_edge[W] = E;
  }
 
  //================================================================
  /// Return pointer to greater or equal-sized edge neighbour
  /// in specified \c direction; also provide info regarding the relative
  /// size and orientation of neighbour:
  /// - The vector translate_s turns the index of the local coordinate
  ///   in the present quadtree into that of the neighbour. If there are no
  ///   rotations then translate_s[i] = i, but if, for example, the neighbour's
  ///   eastern face is adjacent to our northern face translate_s[0] = 1
  ///   and translate_s[1] = 0. Of course, this could be deduced after the
  ///   fact, but it's easier to do it here.
  /// - Let's have a look at the current quadtree's edge in the specified
  ///   direction. The edge will be parallel to one of the two
  ///   local coordinates in the element. This edge has two
  ///   vertices. The vertex at the minimum value of the local
  ///   coordinate (the "lo" vertex) is located
  ///   at the local coordinates (\c s_lo[0], \c s_lo[1]) in the neighbouring
  ///   quadtree.
  /// - ditto with s_hi: The vertex at the maximum value of the local
  ///   coordinate located at the local coordinates (\c s_hi[0],
  ///   \c s_hi[1]) in the neighbouring quadtree.
  /// - We're looking for a neighbour in the specified \c direction. When
  ///   viewed from the neighbouring quadtree, the edge that separates
  ///   the present quadtree from its neighbour is the neighbour's \c edge
  ///   edge. If there's no rotation between the two quadtrees, this is a
  ///   simple reflection: For instance, if we're looking
  ///   for a neighhbour in the \c N [orthern] \c direction, \c edge will
  ///   be \c S [outh]
  /// - \c diff_level <= 0 indicates the difference in refinement levels between
  ///   the two neighbours. If \c diff_level==0, the neighbour has the
  ///   same size as the current quadtree.
  /// - \c in_neighbouring_tree returns true is we have had to flip
  ///      to a different root, even if that root is actually the same
  ///      as it can be in periodic problems.
  //=================================================================
  QuadTree* QuadTree::gteq_edge_neighbour(const int& direction,
                                          Vector<unsigned>& translate_s,
                                          Vector<double>& s_lo,
                                          Vector<double>& s_hi,
                                          int& edge,
                                          int& diff_level,
                                          bool& in_neighbouring_tree) const
  {
    using namespace QuadTreeNames;
 
#ifdef PARANOID
    if ((direction != S) && (direction != E) && (direction != N) &&
        (direction != W))
    {
      std::ostringstream error_stream;
      error_stream << "Direction " << direction << " is not N, S, E, W"
                   << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Initialise in_neighbouring tree to false. It will be set to true
    // during the recursion if we do actually hop over in to the neighbour
    in_neighbouring_tree = false;
 
    // Maximum level to which we're allowed to descend (we only want
    // greater-or-equal-sized neighbours)
    int max_level = Level;
 
    // Current element has the following root:
    QuadTreeRoot* orig_root_pt = dynamic_cast<QuadTreeRoot*>(Root_pt);
 
    // Initialise offset in local coordinate
    double s_diff = 0;
 
    // Initialise difference in level
    diff_level = 0;
 
    // Find neighbour
    QuadTree* return_pt = gteq_edge_neighbour(direction,
                                              s_diff,
                                              diff_level,
                                              in_neighbouring_tree,
                                              max_level,
                                              orig_root_pt);
 
    QuadTree* neighb_pt = return_pt;
 
    // If neighbour exists: What's the direction of the interfacial
    // edge when viewed from within the neighbour element?
    if (neighb_pt != 0)
    {
      // Rotate things around (the orientation of N in the neighbour might be
      // be different from that in the present element)
      // Initialise the direction
      int new_dir = direction;
      // If the neighbour has a different root, then there could be a possible
      // rotation, find it
      if (neighb_pt->root_pt() != Root_pt)
      {
        new_dir = Rotate(orig_root_pt->north_equivalent(direction), direction);
      }
 
      s_lo[0] = S_base(0, Reflect_edge[new_dir]) +
                S_step(0, Reflect_edge[new_dir]) * s_diff;
      s_lo[1] = S_base(1, Reflect_edge[new_dir]) +
                S_step(1, Reflect_edge[new_dir]) * s_diff;
 
      s_hi[0] = S_base(0, Reflect_edge[new_dir]) +
                S_step(0, Reflect_edge[new_dir]) * pow(2.0, diff_level) +
                S_step(0, Reflect_edge[new_dir]) * s_diff;
      s_hi[1] = S_base(1, Reflect_edge[new_dir]) +
                S_step(1, Reflect_edge[new_dir]) * pow(2.0, diff_level) +
                S_step(1, Reflect_edge[new_dir]) * s_diff;
 
      Vector<double> s_lo_new(2), s_hi_new(2);
 
      // What's the direction of the interfacial edge when viewed from within
      // the neighbour element?
      edge = Reflect_edge[new_dir];
 
      // Set up the translation scheme for the local coordinates
      {
        bool swap = false;
        // Do we need to switch the coordinate
        switch (direction)
        {
            // If the direction is north or south,
            // but the neighbour's coordinate s[1] is not constant, we must swap
          case N:
          case S:
            if (s_lo[1] != s_hi[1])
            {
              swap = true;
            }
            break;
            // If the direction is east or west
            // but the neighbour's corodinate s[0] is not constant, we must swap
          case E:
          case W:
            if (s_lo[0] != s_hi[0])
            {
              swap = true;
            }
            break;
          // Catch all totally un-necessary
          default:
            std::ostringstream error_stream;
            error_stream << "Direction " << direction << " is not N, S, E, W"
                         << std::endl;
 
            throw OomphLibError(error_stream.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
        }
 
        // If we must swap, then do so
        if (swap)
        {
          translate_s[0] = 1;
          translate_s[1] = 0;
        }
        // Othewise, it's just a straight translation
        else
        {
          translate_s[0] = 0;
          translate_s[1] = 1;
        }
      }
 
      // Reverse directions?
      s_lo_new[0] = s_lo[0];
      s_lo_new[1] = s_lo[1];
      s_hi_new[0] = s_hi[0];
      s_hi_new[1] = s_hi[1];
 
      if (((edge == N) || (edge == S)) &&
          ((Rotate_angle(direction, new_dir) == 90) ||
           (Rotate_angle(direction, new_dir) == 180)))
      {
        s_lo_new[0] = -s_lo[0];
        s_hi_new[0] = -s_hi[0];
      }
      if (((edge == E) || (edge == W)) &&
          ((Rotate_angle(direction, new_dir) == 270) ||
           (Rotate_angle(direction, new_dir) == 180)))
      {
        s_lo_new[1] = -s_lo[1];
        s_hi_new[1] = -s_hi[1];
      }
 
      s_lo[0] = s_lo_new[0];
      s_lo[1] = s_lo_new[1];
      s_hi[0] = s_hi_new[0];
      s_hi[1] = s_hi_new[1];
    }
    return return_pt;
  }
 
  //================================================================
  /// Find `greater-or-equal-sized edge neighbour' in given direction
  /// (N/E/S/W).
  ///
  /// This is an auxiliary routine which allows neighbour finding in adjacent
  /// quadtrees. Needs to keep track of previous son types and
  /// the maximum level to which search is performed.
  ///
  /// Parameters:
  ///
  /// - direction: N/S/E/W: Direction in which neighbour has to be found.
  /// - s_diff: Offset of lower/left vertex from corresponding vertex in
  ///   neighbour. Note that this is input/output as it needs to be incremented/
  ///   decremented during the recursive calls to this function.
  /// - edge: We're looking for the neighbour across our edge 'direction'
  ///   (N/S/E/W). When viewed from the neighbour, this edge is
  ///   `edge' (N/S/E/W). [If there's no relative rotation between neighbours
  ///    then this is a mere reflection, e.g. direction=N  --> edge=S etc.]
  /// - diff_level <= 0 indicates the difference in quadtree levels
  ///   between the current element and its neighbour.
  /// - max_level is the maximum level to which the neighbour search is
  ///   allowed to proceed. This is again necessary because in a forest,
  ///   the neighbour search isn't based on pure recursion.
  /// - orig_root_pt identifies the root node of the element whose
  ///   neighbour we're really trying to find by all these recursive calls.
  ///
  //=================================================================
  QuadTree* QuadTree::gteq_edge_neighbour(
    const int& direction,
    double& s_diff,
    int& diff_level,
    bool& in_neighbouring_tree,
    int max_level,
    QuadTreeRoot* const& orig_root_pt) const
  {
    using namespace QuadTreeNames;
 
#ifdef PARANOID
    if ((direction != S) && (direction != E) && (direction != N) &&
        (direction != W))
    {
      std::ostringstream error_stream;
      error_stream << "Direction " << direction << " is not N, S, E, W"
                   << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    QuadTree* next_el_pt;
    QuadTree* return_el_pt;
 
    // STEP 1: Find the neighbour's father
    //--------
 
    // Does the element have a father?
    if (Father_pt != 0)
    {
      // If the present quadrant (whose location inside its
      // father element is specified by Son_type) is adjacent to the
      // father's edge in the required direction, then its neighbour has
      // a different father ---> we need to climb up the tree to
      // the father and find his neighbour in the required direction
 
      // Note that this is the cunning recursive part. The returning may not
      // stop until we hit the very top of the tree, when the element does NOT
      // have a father
      if (Is_adjacent(direction, Son_type))
      {
        next_el_pt = dynamic_cast<QuadTree*>(Father_pt)->gteq_edge_neighbour(
          direction,
          s_diff,
          diff_level,
          in_neighbouring_tree,
          max_level,
          orig_root_pt);
      }
      // If the present quadrant is not adjacent to the
      // father's edge in the required direction, then the
      // neighbour has the same father and is obtained
      // by the appropriate reflection inside the father element
      // This will only be called if we have not left the original tree.
      else
      {
        next_el_pt = dynamic_cast<QuadTree*>(Father_pt);
      }
 
      // We're about to ascend one level:
      diff_level -= 1;
 
      // Work out position of lower (or left) corner of present edge
      // in its father element
      s_diff += pow(0.5, -diff_level) * S_direct(direction, Son_type);
 
      // STEP 2:  We have now located the neighbour's father and need to
      // -------  find the appropriate son.
 
      // Buffer cases where the neighbour (and hence its father) lie outside
      // the boundary
      if (next_el_pt != 0)
      {
        // If the father is a leaf then we can't descend to the same
        // level as the present node ---> simply return the father himself
        // as the (greater) neighbour. Same applies if we are about
        // to descend lower than the max_level (in a  neighbouring tree)
        if ((next_el_pt->Son_pt.size() == 0) ||
            (next_el_pt->Level > max_level - 1))
        {
          return_el_pt = next_el_pt;
        }
        // We have located the neighbour's father: The position of the
        // neighbour is obtained by `reflecting' the position of the
        // node itself.
 
        // We know exactly how to reflect, because we know which son type we
        // are and we have the pointer to the neighbours father
        else
        {
          int son_quadrant = Reflect(direction, Son_type);
 
          // If the root of the neighbour's father is not our root, we
          // might need to rotate
          if (orig_root_pt != next_el_pt->Root_pt)
          {
            // Get the north equivalent of the next element
            int my_north =
              dynamic_cast<QuadTreeRoot*>(Root_pt)->north_equivalent(direction);
            son_quadrant = Rotate(my_north, son_quadrant);
          }
 
          // The next element in the tree is the appropriate son of the
          // neighbour's father
          return_el_pt =
            dynamic_cast<QuadTree*>(next_el_pt->Son_pt[son_quadrant]);
 
          // Work out position of lower (or left) corner of present edge
          // in next higher element
          s_diff -= pow(0.5, -diff_level) * S_direct(direction, Son_type);
 
          // We have just descended one level
          diff_level += 1;
        }
      }
      // The neighbour's father lies outside the boundary --> the neighbour
      // itself does too --> return NULL.
      else
      {
        return_el_pt = 0;
      }
    }
    // Element does not have a father --> check if it has a neighbouring
    // tree in the appropriate direction
    else
    {
      // Find  neighbouring root
      if (Root_pt->neighbour_pt(direction) != 0)
      {
        // In this case we have moved to a neighbour, so set the flag
        in_neighbouring_tree = true;
        return_el_pt =
          dynamic_cast<QuadTreeRoot*>(Root_pt->neighbour_pt(direction));
      }
      // No neighbouring tree, so there really is no neighbour --> return NULL
      else
      {
        return_el_pt = 0;
      }
    }
 
    return return_el_pt;
  }
 
  //================================================================
  /// Traverse Tree: Preorder traverse and stick pointers to
  /// neighbouring leaf nodes (only) into Vector
  //=================================================================
  void QuadTree::stick_neighbouring_leaves_into_vector(
    Vector<const QuadTree*>& tree_neighbouring_nodes,
    Vector<Vector<double>>& tree_neighbouring_s_lo,
    Vector<Vector<double>>& tree_neighbouring_s_hi,
    Vector<int>& tree_neighbouring_diff_level,
    const QuadTree* my_neigh_pt,
    const int& direction) const
  {
    // If the tree has sons
    unsigned numsons = Son_pt.size();
    if (numsons > 0)
    {
      // Now do the sons (if they exist)
      for (unsigned i = 0; i < numsons; i++)
      {
        dynamic_cast<QuadTree*>(Son_pt[i])
          ->stick_neighbouring_leaves_into_vector(tree_neighbouring_nodes,
                                                  tree_neighbouring_s_lo,
                                                  tree_neighbouring_s_hi,
                                                  tree_neighbouring_diff_level,
                                                  my_neigh_pt,
                                                  direction);
      }
    }
    else
    {
      // Required data for neighbour-finding routine
      Vector<unsigned> translate_s(2);
      Vector<double> s_lo(2), s_hi(2);
      int edge, diff_level;
      bool in_neighbouring_tree;
      QuadTree* neigh_pt;
 
      // Get neighbouring tree
      neigh_pt = gteq_edge_neighbour(direction,
                                     translate_s,
                                     s_lo,
                                     s_hi,
                                     edge,
                                     diff_level,
                                     in_neighbouring_tree);
 
      // Check if the neighbour is the same as the tree passed in
      // (i.e. Am I a neighbour of the master element's tree?)
      if (neigh_pt == my_neigh_pt)
      {
        // Add the element and the diff_level to passed vectors
        tree_neighbouring_nodes.push_back(this);
        tree_neighbouring_s_lo.push_back(s_lo);
        tree_neighbouring_s_hi.push_back(s_hi);
        tree_neighbouring_diff_level.push_back(diff_level);
      }
    }
  }
 
  //================================================================
  /// Self-test: Check neighbour finding routine. For each element
  /// in the tree and for each vertex, determine the
  /// distance between the vertex and its position in the
  /// neigbour. If the difference is less than
  /// Tree::Max_neighbour_finding_tolerance.
  /// return success (0), otherwise failure (1)
  //=================================================================
  unsigned QuadTree::self_test()
  {
    // Stick pointers to all nodes into Vector and number elements
    // in the process
    Vector<Tree*> all_nodes_pt;
    stick_all_tree_nodes_into_vector(all_nodes_pt);
    long int count = 0;
    unsigned long num_nodes = all_nodes_pt.size();
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      all_nodes_pt[i]->object_pt()->set_number(++count);
    }
 
    // Check neighbours (distance between hanging nodes) -- don't print (keep
    // output streams closed)
    double max_error = 0.0;
    std::ofstream neighbours_file;
    std::ofstream neighbours_txt_file;
    QuadTree::doc_neighbours(
      all_nodes_pt, neighbours_file, neighbours_txt_file, max_error);
 
    if (max_error > Max_neighbour_finding_tolerance)
    {
      oomph_info << "\n \n Failed self_test() for QuadTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 1;
    }
    else
    {
      oomph_info << "\n \n Passed self_test() for QuadTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 0;
    }
  }
 
 
  //================================================================
  /// Constructor for QuadTreeForest:
  ///
  /// Pass:
  ///  - trees_pt[], the Vector of pointers to the constituent trees
  ///    (QuadTreeRoot objects)
  ///
  /// Note that the pointers to the neighbour's of each tree must have
  /// been allocated before the constructor is called, otherwise the
  /// relative rotation scheme will not be constructed correctly.
  //=================================================================
  QuadTreeForest::QuadTreeForest(Vector<TreeRoot*>& trees_pt)
    : TreeForest(trees_pt)
  {
#ifdef LEAK_CHECK
    LeakCheckNames::QuadTreeForest_build += 1;
#endif
 
    // Don't setup neighbours etc. if forest is empty
    if (trees_pt.size() == 0)
    {
      return;
    }
 
    using namespace QuadTreeNames;
 
    // Setup the neighbours
    find_neighbours();
 
    // Construct the rotation scheme, note that all neighbour pointers must
    // be set before the constructor is called
    construct_north_equivalents();
  }
 
 
  //================================================================
  /// Setup the neighbour lookup schemes for all constituent
  /// quadtrees.
  //================================================================
  void QuadTreeForest::find_neighbours()
  {
    using namespace QuadTreeNames;
 
    unsigned numtrees = ntree();
    unsigned n = 0; // to store nnode1d
    if (numtrees > 0)
    {
      n = Trees_pt[0]->object_pt()->nnode_1d();
    }
    else
    {
      throw OomphLibError(
        "Trying to setup the neighbour scheme for an empty forest\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    // Number of vertex nodes: 4
    unsigned n_vertex_node = 4;
 
    // Find potentially connected trees by identifying
    // those whose associated elements share a common vertex node
    std::map<Node*, std::set<unsigned>> tree_assoc_with_vertex_node;
 
    // Loop over all trees
    for (unsigned i = 0; i < numtrees; i++)
    {
      // Loop over the vertex nodes of the associated element
      for (unsigned j = 0; j < n_vertex_node; j++)
      {
        Node* nod_pt = dynamic_cast<QuadElementBase*>(Trees_pt[i]->object_pt())
                         ->vertex_node_pt(j);
        tree_assoc_with_vertex_node[nod_pt].insert(i);
      }
    }
 
 
    // For each tree we store a set of potentially neighbouring trees
    // i.e. trees that share at least one node
    Vector<std::set<unsigned>> potentially_neighb_tree(numtrees);
 
    // Loop over vertex nodes
    for (std::map<Node*, std::set<unsigned>>::iterator it =
           tree_assoc_with_vertex_node.begin();
         it != tree_assoc_with_vertex_node.end();
         it++)
    {
      // Loop over connected elements twice
      for (std::set<unsigned>::iterator it_el1 = it->second.begin();
           it_el1 != it->second.end();
           it_el1++)
      {
        unsigned i = (*it_el1);
        for (std::set<unsigned>::iterator it_el2 = it->second.begin();
             it_el2 != it->second.end();
             it_el2++)
        {
          unsigned j = (*it_el2);
          // These two elements are potentially connected
          if (i != j)
          {
            potentially_neighb_tree[i].insert(j);
          }
        }
      }
    }
 
 
    // Loop over all trees
    for (unsigned i = 0; i < numtrees; i++)
    {
      // Loop over their potential neighbours
      for (std::set<unsigned>::iterator it = potentially_neighb_tree[i].begin();
           it != potentially_neighb_tree[i].end();
           it++)
      {
        unsigned j = (*it);
 
        // is it the Northern neighbour ?
        bool is_N_neighbour =
          ((Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1) &&
           (Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
 
 
        // is it the Southern neighbour ?
        bool is_S_neighbour =
          ((Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
           (Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1));
 
 
        // is it the Eastern neighbour ?
        bool is_E_neighbour =
          ((Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1) &&
           (Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
 
        // is it the Western neighbour ?
        bool is_W_neighbour =
          ((Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
           (Trees_pt[j]->object_pt()->get_node_number(
              Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1));
 
 
        if (is_N_neighbour) Trees_pt[i]->neighbour_pt(N) = Trees_pt[j];
        if (is_S_neighbour) Trees_pt[i]->neighbour_pt(S) = Trees_pt[j];
        if (is_E_neighbour) Trees_pt[i]->neighbour_pt(E) = Trees_pt[j];
        if (is_W_neighbour) Trees_pt[i]->neighbour_pt(W) = Trees_pt[j];
      }
    }
 
 
    // Old hacky version with horrendous scaling
    if (false)
    {
      // Loop over all trees
      for (unsigned i = 0; i < numtrees; i++)
      {
        // Loop over all the other elements
        for (unsigned j = 0; j < numtrees; j++)
        {
          if (j != i) // make sure we are not looking at the element itself
          {
            // is it the Northern neighbour ?
            bool is_N_neighbour =
              ((Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1) &&
               (Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
 
 
            // is it the Southern neighbour ?
            bool is_S_neighbour =
              ((Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
               (Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1));
 
 
            // is it the Eastern neighbour ?
            bool is_E_neighbour =
              ((Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1) &&
               (Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
 
            // is it the Western neighbour ?
            bool is_W_neighbour =
              ((Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
               (Trees_pt[j]->object_pt()->get_node_number(
                  Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1));
 
 
            if (is_N_neighbour) Trees_pt[i]->neighbour_pt(N) = Trees_pt[j];
            if (is_S_neighbour) Trees_pt[i]->neighbour_pt(S) = Trees_pt[j];
            if (is_E_neighbour) Trees_pt[i]->neighbour_pt(E) = Trees_pt[j];
            if (is_W_neighbour) Trees_pt[i]->neighbour_pt(W) = Trees_pt[j];
          }
        }
      }
    }
  }
 
 
  //================================================================
  /// Construct the rotation scheme for the quadtree forest.
  /// Note that all pointers to neighbours must have been allocated
  /// for this to work.
  //================================================================
  void QuadTreeForest::construct_north_equivalents()
  {
    using namespace QuadTreeNames;
 
    unsigned numtrees = ntree();
    // Loop over all the trees
    for (unsigned i = 0; i < numtrees; i++)
    {
      // Find the pointer to the northern neighbour
      QuadTreeRoot* neigh_pt = quad_neigh_pt(i, N);
      // If there is a neighbour
      if (neigh_pt != 0)
      {
        // Find the direction of the present tree, as viewed from the neighbour
        int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
 
        // Set up the rotation scheme
        switch (direction)
        {
            // If N neighbour has this tree on S, north equivalent is N
          case S:
            quadtree_pt(i)->north_equivalent(N) = N;
            break;
            // If N neighbour has this tree on W, north equivalent is E
          case W:
            quadtree_pt(i)->north_equivalent(N) = E;
            break;
            // If N neighbour has this tree on N, north equivalent is S
          case N:
            quadtree_pt(i)->north_equivalent(N) = S;
            break;
            // If N neighbour has this tree on E, north equivalent is W
          case E:
            quadtree_pt(i)->north_equivalent(N) = W;
            break;
            // If N neighbour does not have pointer to this tree, die
          default:
            std::ostringstream error_stream;
            error_stream
              << "Tree " << i
              << "'s Northern neighbour has no neighbour pointer to Tree " << i
              << std::endl;
 
            throw OomphLibError(error_stream.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Find the pointer to the eastern neighbour
      neigh_pt = quad_neigh_pt(i, E);
      // If there is a neighbour
      if (neigh_pt != 0)
      {
        // Find the direction of the present tree, as viewed from the neighbour
        int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
 
        // Set up the rotation scheme
        switch (direction)
        {
            // If E neighbour has this tree on W, north equivalent is N
          case W:
            quadtree_pt(i)->north_equivalent(E) = N;
            break;
            // If E neighbour has this tree on N, north equivalent is E
          case N:
            quadtree_pt(i)->north_equivalent(E) = E;
            break;
            // If E neighbour has this tree on E, north equivalent is S
          case E:
            quadtree_pt(i)->north_equivalent(E) = S;
            break;
            // If E neighbour has this tree on S, north equivalent is W
          case S:
            quadtree_pt(i)->north_equivalent(E) = W;
            break;
            // If E neighbour does not have pointer to this tree, die
          default:
            std::ostringstream error_stream;
            error_stream
              << "Tree " << i
              << "'s Eastern neighbour has no neighbour pointer to Tree " << i
              << std::endl;
 
            throw OomphLibError(error_stream.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Find the pointer to the southern neighbour
      neigh_pt = quad_neigh_pt(i, S);
      // If there is a neighbour
      if (neigh_pt != 0)
      {
        // Find the direction of the present tree, as viewed from the neighbour
        int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
 
        // Set up the rotation scheme
        switch (direction)
        {
            // If S neighbour has this tree on N, north equivalent is N
          case N:
            quadtree_pt(i)->north_equivalent(S) = N;
            break;
            // If S neighbour has this tree on E, north equivalent is E
          case E:
            quadtree_pt(i)->north_equivalent(S) = E;
            break;
            // If S neighbour has this tree on S, north equivalent is S
          case S:
            quadtree_pt(i)->north_equivalent(S) = S;
            break;
            // If S neighbour has this tree on W, north equivalent is W
          case W:
            quadtree_pt(i)->north_equivalent(S) = W;
            break;
            // If S neighbour does not have pointer to this tree, die
          default:
            std::ostringstream error_stream;
            error_stream
              << "Tree " << i
              << "'s Southern neighbour has no neighbour pointer to Tree " << i
              << std::endl;
 
            throw OomphLibError(error_stream.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Find the pointer to the western neighbour
      neigh_pt = quad_neigh_pt(i, W);
      // If there is a neighbour
      if (neigh_pt != 0)
      {
        // Find the direction of the present tree, as viewed from the neighbour
        int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
 
        // Set up the rotation scheme
        switch (direction)
        {
            // If W neighbour has this tree on E, north equivalent is N
          case E:
            quadtree_pt(i)->north_equivalent(W) = N;
            break;
            // If W neighbour has this tree on S, north equivalent is E
          case S:
            quadtree_pt(i)->north_equivalent(W) = E;
            break;
            // If W neighbour has this tree on W, north equivalent is S
          case W:
            quadtree_pt(i)->north_equivalent(W) = S;
            break;
            // If W neighbour has this tree on N, north equivalent is W
          case N:
            quadtree_pt(i)->north_equivalent(W) = W;
            break;
            // If W neighbour does not have pointer to this tree, die
          default:
            std::ostringstream error_stream;
            error_stream
              << "Tree " << i
              << "'s Western neighbour has no neighbour pointer to Tree " << i
              << std::endl;
 
            throw OomphLibError(error_stream.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
        }
      }
    }
  }
 
  //================================================================
  /// Document and check all the neighbours in all the nodes in
  /// the forest
  //================================================================
  void QuadTreeForest::check_all_neighbours(DocInfo& doc_info)
  {
    // Create Vector of elements
    Vector<Tree*> all_tree_nodes_pt;
    this->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
 
    // Create storage for information files
    std::ofstream neigh_file;
    std::ofstream neigh_txt_file;
 
    // If we are documenting the results, then open the files
    if (doc_info.is_doc_enabled())
    {
      std::ostringstream fullname;
      fullname << doc_info.directory() << "/neighbours" << doc_info.number()
               << ".dat";
      oomph_info << "opened " << fullname.str() << " to doc neighbours"
                 << std::endl;
      neigh_file.open(fullname.str().c_str());
      fullname.str("");
      fullname << doc_info.directory() << "/neighbours" << doc_info.number()
               << ".txt";
      oomph_info << "opened " << fullname.str() << " to doc neighbours"
                 << std::endl;
      neigh_txt_file.open(fullname.str().c_str());
    }
 
    // Call the standard documentation function
    double max_error = 0.0;
    QuadTree::doc_neighbours(
      all_tree_nodes_pt, neigh_file, neigh_txt_file, max_error);
 
    // If the error is too large complain
    if (max_error > Tree::max_neighbour_finding_tolerance())
    {
      std::ostringstream error_stream;
      error_stream << "Max. error in quadtree neighbour finding: " << max_error
                   << " is too big" << std::endl;
      error_stream
        << "i.e. bigger than Tree::max_neighbour_finding_tolerance()="
        << Tree::max_neighbour_finding_tolerance() << std::endl;
 
      // Close the files if they were opened
      if (doc_info.is_doc_enabled())
      {
        neigh_file.close();
        neigh_txt_file.close();
      }
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    else
    {
      oomph_info << "Max. error in quadtree neighbour finding: " << max_error
                 << " is OK" << std::endl;
      oomph_info
        << "i.e. less than QuadTree::max_neighbour_finding_tolerance()="
        << QuadTree::max_neighbour_finding_tolerance() << std::endl;
    }
 
    // Close the files if they were opened
    if (doc_info.is_doc_enabled())
    {
      neigh_file.close();
      neigh_txt_file.close();
    }
  }
 
  //================================================================
  /// Open output files that will stored any hanging nodes that are
  /// created in the mesh refinement process.
  /// ===============================================================
  void QuadTreeForest::open_hanging_node_files(
    DocInfo& doc_info, Vector<std::ofstream*>& output_stream)
  {
    // In 2D, there will be four output files
    for (unsigned i = 0; i < 4; i++)
    {
      output_stream.push_back(new std::ofstream);
    }
 
    // If we are documenting the output, open the files
    if (doc_info.is_doc_enabled())
    {
      std::ostringstream fullname;
      fullname << doc_info.directory() << "/hang_nodes_s" << doc_info.number()
               << ".dat";
      output_stream[0]->open(fullname.str().c_str());
      fullname.str("");
      fullname << doc_info.directory() << "/hang_nodes_n" << doc_info.number()
               << ".dat";
      output_stream[1]->open(fullname.str().c_str());
      fullname.str("");
      fullname << doc_info.directory() << "/hang_nodes_w" << doc_info.number()
               << ".dat";
      output_stream[2]->open(fullname.str().c_str());
      fullname.str("");
      fullname << doc_info.directory() << "/hang_nodes_e" << doc_info.number()
               << ".dat";
      output_stream[3]->open(fullname.str().c_str());
    }
  }
 
  //================================================================
  /// Self test:  Check neighbour finding routine. For each element
  /// in the tree and for each vertex, determine the
  /// distance between the vertex and its position in the
  /// neigbour. If the difference is less than
  /// Tree::Max_neighbour_finding_tolerance.
  /// return success (0), otherwise failure (1)
  //=================================================================
  unsigned QuadTreeForest::self_test()
  {
    // Stick pointers to all nodes into Vector and number elements
    // in the process
    Vector<Tree*> all_forest_nodes_pt;
    stick_all_tree_nodes_into_vector(all_forest_nodes_pt);
    long int count = 0;
    unsigned long num_nodes = all_forest_nodes_pt.size();
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      all_forest_nodes_pt[i]->object_pt()->set_number(++count);
    }
 
    // Check neighbours (distance between hanging nodes) -- don't print
    // (keep output streams closed)
    double max_error = 0.0;
    std::ofstream neighbours_file;
    std::ofstream neighbours_txt_file;
    QuadTree::doc_neighbours(
      all_forest_nodes_pt, neighbours_file, neighbours_txt_file, max_error);
    if (max_error > QuadTree::max_neighbour_finding_tolerance())
    {
      oomph_info << "\n \n Failed self_test() for QuadTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 1;
    }
    else
    {
      oomph_info << "\n \n Passed self_test() for QuadTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 0;
    }
  }
 
  //=================================================================
  /// Doc/check all neighbours of quadtree ("nodes") contained in the
  /// Vector forest_node_pt. Output into neighbours_file which can
  /// be viewed from tecplot with QuadTreeNeighbours.mcr
  /// Neighbour info and errors are displayed on
  /// neighbours_txt_file.  Finally, compute the max. error between
  /// vertices when viewed from neighhbouring element.
  /// Output is suppressed if the output streams are closed.
  //=================================================================
  void QuadTree::doc_neighbours(Vector<Tree*> forest_nodes_pt,
                                std::ofstream& neighbours_file,
                                std::ofstream& neighbours_txt_file,
                                double& max_error)
  {
    using namespace QuadTreeNames;
 
    int diff_level;
    double s_diff;
    bool in_neighbouring_tree;
    int edge = OMEGA;
 
    Vector<double> s(2);
    Vector<double> x(2);
    Vector<int> prev_son_type;
 
    Vector<unsigned> translate_s(2);
    Vector<double> s_lo(2);
    Vector<double> s_hi(2);
 
    Vector<double> x_small(2);
    Vector<double> x_large(2);
 
    // Initialise error in vertex positions
    max_error = 0.0;
 
 
    // Loop over all elements to assign numbers for plotting
    // -----------------------------------------------------
    unsigned long num_nodes = forest_nodes_pt.size();
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      // Set number
      forest_nodes_pt[i]->object_pt()->set_number(i);
    }
 
    // Loop over all elements for checks
    // ---------------------------------
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      // Doc the element itself
      QuadTree* el_pt = dynamic_cast<QuadTree*>(forest_nodes_pt[i]);
 
      // If the object is incomplete complain
      if (el_pt->object_pt()->nodes_built())
      {
        // Print it
        if (neighbours_file.is_open())
        {
          dynamic_cast<RefineableQElement<2>*>(el_pt->object_pt())
            ->output_corners(neighbours_file, "BLACK");
        }
 
        // Loop over directions to find neighbours
        // ----------------------------------------
        for (int direction = N; direction <= W; direction++)
        {
          // Initialise difference in levels and coordinate offset
          diff_level = 0;
          s_diff = 0.0;
 
          // Find greater-or-equal-sized neighbour...
          QuadTree* neighb_pt =
            el_pt->gteq_edge_neighbour(direction,
                                       translate_s,
                                       s_lo,
                                       s_hi,
                                       edge,
                                       diff_level,
                                       in_neighbouring_tree);
 
          // If neighbour exist and nodes are created: Doc it
          if ((neighb_pt != 0) && (neighb_pt->object_pt()->nodes_built()))
          {
            // Doc neighbour stats
            if (neighbours_txt_file.is_open())
            {
              neighbours_txt_file
                << Direct_string[direction] << " neighbour of "
                << el_pt->object_pt()->number() << " is "
                << neighb_pt->object_pt()->number() << " diff_level "
                << diff_level << " s_diff " << s_diff
                << " inside neighbour the edge is " << Direct_string[edge]
                << std::endl
                << std::endl;
            }
 
            // Plot neighbour in the appropriate Colour
            if (neighbours_file.is_open())
            {
              dynamic_cast<RefineableQElement<2>*>(neighb_pt->object_pt())
                ->output_corners(neighbours_file, Colour[direction]);
            }
 
            // Check that local
            // coordinates in the larger element (obtained via s_diff)
            // lead to the same spatial point as the node vertices
            // in the current element
            {
              if (neighbours_file.is_open())
              {
                neighbours_file << "ZONE I=2 \n";
              }
 
              // (Lower)/(left) vertex:
              //-----------------------
 
              // Get coordinates in large (neighbour) element
              s[0] = s_lo[0];
              s[1] = s_lo[1];
              neighb_pt->object_pt()->get_x(s, x_large);
 
              // Get coordinates in small element
              Vector<double> s(2);
              s[0] = S_base(0, direction);
              s[1] = S_base(1, direction);
              el_pt->object_pt()->get_x(s, x_small);
 
 
              // Need to exclude periodic nodes from this check
              // There can only be periodic nodes if we have moved into the
              // neighbour
              bool is_periodic = false;
              if (in_neighbouring_tree)
              {
                // is the node periodic
                is_periodic =
                  el_pt->root_pt()->is_neighbour_periodic(direction);
              }
 
              double error = 0.0;
              // Only bother to calculate the error if the node is NOT periodic
              if (is_periodic == false)
              {
                for (int i = 0; i < 2; i++)
                {
                  error += pow(x_small[i] - x_large[i], 2);
                }
              }
 
              // Take the root of the square error
              error = sqrt(error);
              if (neighbours_txt_file.is_open())
              {
                neighbours_txt_file << "Error (1) " << error << std::endl;
              }
 
              if (std::fabs(error) > max_error)
              {
                max_error = std::fabs(error);
              }
 
              if (neighbours_file.is_open())
              {
                neighbours_file << x_large[0] << " " << x_large[1] << "  0 \n";
              }
 
              // (Upper)/(right) vertex:
              //------------------------
 
              // Get coordinates in large (neighbour) element
              s[0] = s_hi[0];
              s[1] = s_hi[1];
              neighb_pt->object_pt()->get_x(s, x_large);
 
              // Get coordinates in small element
              s[0] = S_base(0, direction) + S_step(0, direction);
              s[1] = S_base(1, direction) + S_step(1, direction);
              el_pt->object_pt()->get_x(s, x_small);
 
              error = 0.0;
              // Only do this if we are NOT periodic
              if (is_periodic == false)
              {
                for (int i = 0; i < 2; i++)
                {
                  error += pow(x_small[i] - x_large[i], 2);
                }
              }
              // Take the root of the square error
              error = sqrt(error);
 
              // error=
              // sqrt(pow(x_small[0]-x_large[0],2)+pow(x_small[1]-x_large[1],2));
              if (neighbours_txt_file.is_open())
              {
                neighbours_txt_file << "Error (2) " << error << std::endl;
              }
 
              if (std::fabs(error) > max_error)
              {
                max_error = std::fabs(error);
              }
 
              if (neighbours_file.is_open())
              {
                neighbours_file << x_large[0] << " " << x_large[1] << "  0 \n";
              }
            }
            //        else
            //         {
            //          // No neighbour: write dummy zone so tecplot can find
            //          four
            //          // neighbours for every element
            //          if (neighbours_file.is_open())
            //            {
            //             neighbours_file << "ZONE I=2 \n";
            //             neighbours_file << "-0.05 -0.05   0 \n";
            //             neighbours_file << "-0.05 -0.05   0 \n";
            //            }
            //         }
          }
          // If neighbour does not exist: Insert blank zones into file
          // so that tecplot can find four neighbours for every element
          else
          {
            if (neighbours_file.is_open())
            {
              neighbours_file << "ZONE \n 0.00 0.00  0 \n";
              neighbours_file << "ZONE I=2 \n";
              neighbours_file << "-0.05 -0.05  0 \n";
              neighbours_file << "-0.05 -0.05  0 \n";
            }
          }
        }
      } // End of case when element can be documented
    }
  }
 
} // namespace oomph