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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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
  }
}

}