octree.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//====================================================================
#include<algorithm>

#include "octree.h"
#include "refineable_brick_element.h"


namespace oomph
{

//========================================================================
/// Bool indicating that static member data has been setup
//========================================================================
bool OcTree::Static_data_has_been_setup=false;



// Public static member data:
//---------------------------

//====================================================================
/// Translate (enumerated) directions into strings 
//====================================================================
Vector<std::string> OcTree::Direct_string;

//====================================================================
/// Get opposite face, e.g. Reflect_face[L]=R 
//====================================================================
Vector<int> OcTree::Reflect_face;

//====================================================================
/// Get opposite edge, e.g. Reflect_edge[DB]=UF 
//====================================================================
Vector<int> OcTree::Reflect_edge;

//====================================================================
/// Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF 
//====================================================================
Vector<int> OcTree::Reflect_vertex;

//====================================================================
/// Map of vectors containing the two vertices for each edge
//====================================================================
Vector<Vector<int> > OcTree::Vertex_at_end_of_edge;

//====================================================================
/// Map storing the information to translate a vector of directions 
/// (in the three coordinate directions) into a direction
//====================================================================
std::map<Vector<int>, int > OcTree::Vector_to_direction;

//====================================================================
/// Vector storing the information to translate a direction into a
/// vector of directions (in the three coordinate directions)
//====================================================================
Vector<Vector<int> > OcTree::Direction_to_vector;

//====================================================================
/// \short Storage for the up/right-equivalents corresponding to two
/// pairs of vertices along an element edge:
/// - The first pair contains 
///   -# the vertex in the reference element
///   -# the corresponding vertex in the edge neighbour (i.e. the
///      vertex in the edge neighbour that is located at the same
///      position as that first vertex). 
///   .
/// - The second pair contains
///    -# the vertex at the other end of the edge in the reference element
///    -# the corresponding vertex in the edge neighbour.
///    .
/// .
/// These two pairs completely define the relative rotation 
/// between the reference element and its edge neighbour. The map
/// returns a pair which contains the up_equivalent and the
/// right_equivalent of the edge neighbour, i.e. it tells us
/// which direction in the edge neighbour coincides with the
/// up (or right) direction in the reference element. 
//====================================================================
std::map<std::pair<std::pair<int,int>,std::pair<int,int> >,std::pair<int,int> >
OcTree::Up_and_right_equivalent_for_pairs_of_vertices;





// Protected static member data:
//------------------------------


//====================================================================
 /// Entry in rotation matrix: cos(i*90)
//====================================================================
Vector<int> OcTree::Cosi;


//====================================================================
 /// Entry in rotation matrix: sin(i*90)
//====================================================================
Vector<int> OcTree::Sini; 

//====================================================================
/// Array of direction/octant adjacency scheme: 
/// Is_adjacent(direction,octant): Is face/edge \c direction
/// adjacent to octant \c octant ? Table in Samet's book.
//==================================================================== 
DenseMatrix<bool> OcTree::Is_adjacent;

//====================================================================
/// Reflection scheme: Reflect(direction,octant): Get mirror 
/// of octant/edge in specified direction. E.g. Reflect(LDF,L)=RDF
//====================================================================
DenseMatrix<int> OcTree::Reflect;

//====================================================================
/// Determine common face of edges or octants.
/// Slightly bizarre lookup scheme from Samet's book.
//====================================================================
DenseMatrix<int> OcTree::Common_face;

//====================================================================
/// Colours for neighbours in various directions
//====================================================================
Vector<std::string> OcTree::Colour;

//====================================================================
/// s_base(i,direction):  Initial value for coordinate s[i] on 
/// the face indicated by direction (L/R/U/D/F/B) 
//====================================================================
DenseMatrix<double> OcTree::S_base;

//====================================================================
/// Each face of the RefineableQElement<3> that is represented 
/// by the octree is parametrised by two (of the three) 
/// local coordinates that parametrise the entire 3D element. E.g. 
/// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face 
/// is parametrised by (s[0],s[2]); etc. We always identify the 
/// in-face coordinate with the lower (3D) index with the subscript 
/// \c _lo and the one with the larger (3D) index with the subscript \c _hi.
///  Here we set up the translation scheme between the 2D in-face
/// coordinates (s_lo,s_hi) and the corresponding 3D coordinates:
/// If we're located on face \c face [L/R/F/B/U/D], then
/// an increase in s_lo from -1 to +1 corresponds to a change
/// of \c S_steplo(i,face) in the 3D coordinate \c s[i].  S_steplo(i,direction)
//====================================================================
DenseMatrix<double> OcTree::S_steplo; 


//====================================================================
/// If we're located on face \c face [L/R/F/B/U/D], then
/// an increase in s_hi from -1 to +1 corresponds to a change
/// of \c S_stephi(i,face) in the 3D coordinate \ s[i]. 
/// [Read the discussion of \c S_steplo for an explanation of
/// the subscripts \c _hi and \c _lo.]
//====================================================================
DenseMatrix<double> OcTree::S_stephi; 

//====================================================================
///  Relative to the left/down/back vertex in any (father) octree, the 
/// corresponding vertex in the son specified by \c son_octant has an offset.
/// If we project the son_octant's left/down/back vertex onto the
/// father's face \c face, it is located at the in-face coordinate 
///  \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
///  \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
//====================================================================
DenseMatrix<double> OcTree::S_directlo;

//====================================================================
/// Relative to the left/down/back vertex in any (father) octree, the 
/// corresponding vertex in the son specified by \c son_octant has an offset.
/// If we project the son_octant's left/down/back vertex onto the
/// father's face \c face, it is located at the in-face coordinate 
/// \c s_hi = h/2 \c S_directlhi(face,son_octant). [See discussion of
/// \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
//====================================================================
DenseMatrix<double> OcTree::S_directhi;

//====================================================================
/// S_base_edge(i,edge):  Initial value for coordinate s[i] on 
/// the specified edge (LF/RF/...).
//====================================================================
DenseMatrix<double> OcTree::S_base_edge;

//====================================================================
///  Each edge of the RefineableQElement<3> that is represented 
 /// by the octree  is parametrised by one (of the three) 
 /// local coordinates that parametrise the entire 3D element.
 /// If we're located on edge \c edge [DB,UB,...], then
 /// an increase in s from -1 to +1 corresponds to a change
 /// of \c s_step_edge(i,edge) in the 3D coordinates \c s[i]. 
//====================================================================
DenseMatrix<double> OcTree::S_step_edge; 

//====================================================================
///  Relative to the left/down/back vertex in any (father) octree, the 
/// corresponding vertex in the son specified by \c son_octant has an offset.
/// If we project the son_octant's left/down/back vertex onto the
/// father's edge \c edge, it is located at the in-face coordinate 
///  \c s_lo = h/2 \c S_direct_edge(edge,son_octant). 
//====================================================================
DenseMatrix<double> OcTree::S_direct_edge;





// End static data
//----------------



//===================================================================
/// This function is used to translate the position of a vertex node
/// (given by his local number n into a vector giving the position of
/// this node in the local coordinates system.
/// It also needs the value of nnode1d to work.
//===================================================================
Vector<int> OcTree::vertex_node_to_vector(const unsigned& n, 
                                          const unsigned& nnode1d)
{
#ifdef PARANOID
 if ((n!=0)&&
     (n!=nnode1d-1)&&
     (n!=(nnode1d-1)*nnode1d)&&
     (n!=nnode1d*nnode1d-1)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+0)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d-1)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+(nnode1d-1)*nnode1d)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d*nnode1d-1))
  {
   std::ostringstream error_stream;
   error_stream << "Node " << n 
                << " is not a vertex node in a brick element with " 
                << nnode1d << " nodes along each edge!";
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
 int a,b,c;
 Vector<int> result_vect(3);
 a=n / (nnode1d*nnode1d);
 b=(n-a*nnode1d*nnode1d)/nnode1d;
 c=n-a*nnode1d*nnode1d-b*nnode1d;

 result_vect[0]=2*c/(nnode1d-1)-1;
 result_vect[1]=2*b/(nnode1d-1)-1;
 result_vect[2]=2*a/(nnode1d-1)-1;
 
 return result_vect;
}





//==================================================================
/// Return the local node number of given vertex
/// in an element with nnode1d nodes in each coordinate direction.
//==================================================================
unsigned OcTree::vertex_to_node_number(const int& vertex, 
                                       const unsigned& nnode1d)
{
 using namespace OcTreeNames;

#ifdef PARANOID
    if ((vertex!=LDB)&&(vertex!=RDB)&&
        (vertex!=LUB)&&(vertex!=RUB)&&
        (vertex!=LDF)&&(vertex!=RDF)&&
        (vertex!=LUF)&&(vertex!=RUF))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong vertex: " << Direct_string[vertex] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

 switch (vertex)
  {

  case LDB:
   return 0;
   break;

  case RDB:
   return nnode1d-1;
   break;


  case LUB:
   return nnode1d*(nnode1d-1);
   break;

  case RUB:
   return nnode1d*nnode1d-1;
   break;


  case LDF:
   return nnode1d*nnode1d*(nnode1d-1);
   break;


  case RDF:
   return (nnode1d*nnode1d+1)*(nnode1d-1);
   break;

  case LUF:
   return nnode1d*nnode1d*nnode1d-nnode1d;
   break;

  case RUF:
   return nnode1d*nnode1d*nnode1d-1;
   break;
   
  default:
   
   std::ostringstream error_stream;
   error_stream << "Never get here. vertex: " << Direct_string[vertex] 
                << std::endl;
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
   break;
  }

        
        
 std::ostringstream error_stream;
 error_stream << "Never get here. vertex: " << Direct_string[vertex] 
              << std::endl;
 throw OomphLibError(error_stream.str(),
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);


}


//==================================================================
/// \short Return the vertex of local (vertex) node n
/// in an element with nnode1d nodes in each coordinate direction.
//==================================================================
 int OcTree::node_number_to_vertex(const unsigned& n, 
                                   const unsigned& nnode1d)
{
 using namespace OcTreeNames;
 
#ifdef PARANOID
 if ((n!=0)&&
     (n!=nnode1d-1)&&
     (n!=(nnode1d-1)*nnode1d)&&
     (n!=nnode1d*nnode1d-1)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+0)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d-1)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+(nnode1d-1)*nnode1d)&&
     (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d*nnode1d-1))
  {
   std::ostringstream error_stream;
   error_stream << "Node " << n 
                << " is not a vertex node in a brick element with " 
                << nnode1d << " nodes along each edge!";
   
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
 if (n==0) {return LDB;}
 else if (n==nnode1d-1) {return RDB;}
 else if (n==nnode1d*(nnode1d-1)) {return LUB;}
 else if (n==nnode1d*nnode1d-1) {return RUB;}
 else if (n==nnode1d*nnode1d*(nnode1d-1)) {return LDF;}
 else if (n==(nnode1d*nnode1d+1)*(nnode1d-1)) {return RDF;}
 else if (n==nnode1d*nnode1d*nnode1d-nnode1d) {return LUF;}
 else if (n==nnode1d*nnode1d*nnode1d-1) {return RUF;}
 else
  {
   std::ostringstream error_stream;
   error_stream << "Never get here. local node number: " << n
                << " is not a vertex node in a brick element with " 
                << nnode1d << " nodes along each edge!"
                << std::endl;
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
} 




//==================================================================
/// This function takes as argument two node numbers of two nodes
/// delimiting an edge, and one face of this edge and returns the 
/// other face that is sharing this edge. The node numbers given to
/// this function MUST be vertices nodes to work.
/// it also need the value of nnode1d to work.
/// (\c face is a direction in the set U,D,F,B,L,R).
//==================================================================
int OcTree::get_the_other_face(const unsigned& n1, 
                               const unsigned& n2, 
                               const unsigned& nnode1d,
                               const int& face)
{
 Vector<int> vect_node1(3);
 Vector<int> vect_node2(3);
 Vector<int> vect_face(3);
 Vector<int> vect_other_face(3);
 
 // translate the nodes to vectors
 vect_node1=vertex_node_to_vector(n1,nnode1d);
 vect_node2=vertex_node_to_vector(n2,nnode1d);
 //translate the face to a vector
 vect_face=Direction_to_vector[face];

 // compute the vector to the other face -- magic, courtesy of Renaud
 // Schleck!
 for(unsigned i=0;i<3;i++)
  {
   vect_other_face[i]=(vect_node1[i]+vect_node2[i])/2-vect_face[i];
   
#ifdef PARANOID
   if((vect_other_face[i]!=1)&&(vect_other_face[i]!=-1)&&
      (vect_other_face[i]!=0))
    {
     throw OomphLibError("The nodes given are not vertices",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
  }
  
 //return the corresponding face
 return Vector_to_direction[vect_other_face];
}



//====================================================================
/// Build the rotation matrix for a rotation around the axis \c axis of
/// an angle \c angle*90
//====================================================================
void OcTree::construct_rotation_matrix(int& axis, int& angle,
                                       DenseMatrix<int>& mat)
{
  using namespace OcTreeNames;
  int a=0,b=0,c=0,i,j;

  switch(axis)
    {
    case R : a=1; b=2; c=0; break;
    case U : a=2; b=0; c=1; break;
    case F : a=0; b=1; c=2; break;
    default : 
     throw OomphLibError("Bad Axis",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
  for(i=0;i<3;i++)
    {
      for(j=0;j<3;j++)
        {
          mat(i,j)=0;
        }
    }
  mat(a,a) = Cosi[angle];
  mat(b,b) = Cosi[angle];
  mat(a,b) = -Sini[angle];
  mat(b,a) = Sini[angle];
  mat(c,c) = 1;
} 

//===================================================================
/// Helper: Performs the operation Mat3=Mat1*Mat2
//===================================================================
void OcTree::mult_mat_mat(const DenseMatrix<int>& mat1,
                          const DenseMatrix<int>& mat2,
                          DenseMatrix<int>& mat3) 
{
  int Sum,i,j,k;
  for(i=0;i<3;i++)
    {
      for(j=0;j<3;j++)
        {
          Sum=0;
          for(k=0;k<3;k++)
            {
              Sum+=mat1(i,k)*mat2(k,j);
            }
          mat3(i,j)=Sum;
         }
    }
}



//===================================================================
/// Helper: Performs the operation Vect2=Mat*Vect1
//===================================================================
void OcTree::mult_mat_vect(const DenseMatrix<int>& mat,
                           const Vector<int>& vect1, 
                           Vector<int>& vect2) 
{
  int i,k,sum;
  for(i=0;i<3;i++)
    {
      sum=0;
      for(k=0;k<3;k++)
        {
          sum+=mat(i,k)*vect1[k];
        }
      vect2[i]=sum;
    }
}




//===================================================================
/// A rotation is defined by the newUp and newRight
/// directions; so if Up becomes newUp and Right becomes newRight
/// then dir becomes rotate(newUp,newRight,dir);
//===================================================================
int OcTree::rotate(const int& new_up, const int& new_right,
                   const int& dir) 
{
 using namespace OcTreeNames;

#ifdef PARANOID
    if ((new_up!=L)&&(new_up!=R)&&
        (new_up!=F)&&(new_up!=B)&&
        (new_up!=U)&&(new_up!=D))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong new_up: " << Direct_string[new_up] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
    if ((new_right!=L)&&(new_right!=R)&&
        (new_right!=F)&&(new_right!=B)&&
        (new_right!=U)&&(new_right!=D))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong new_right: " << Direct_string[new_right] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

 Vector<int> vect_dir(3);
 Vector<int> vect_new_dir(3);
 // translate the direction to a vector
 vect_dir=Direction_to_vector[dir];
 
 //rotate it
 vect_new_dir=rotate(new_up,new_right,vect_dir);
 
 // translate the vector back to a direction
 return Vector_to_direction[vect_new_dir];
}


//==================================================================
/// This function rotates a vector according to a rotation of
/// the axes that changes up to new_up and right to new_right.
//==================================================================
Vector<int> OcTree::rotate(const int& new_up,const int& new_right,
                           const Vector<int>& dir) 
{
  using namespace OcTreeNames;
  
#ifdef PARANOID
    if ((new_up!=L)&&(new_up!=R)&&
        (new_up!=F)&&(new_up!=B)&&
        (new_up!=U)&&(new_up!=D))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong new_up: " << Direct_string[new_up] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
    if ((new_right!=L)&&(new_right!=R)&&
        (new_right!=F)&&(new_right!=B)&&
        (new_right!=U)&&(new_right!=D))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong new_right: " << Direct_string[new_right] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif


  // Every possible rotation of an element can be produced by the 
  // compostition of two (or one) rotations about one of the main
  // axes of the element (R, U, F). So for each (newRight,newUp)
  // we define the (angle1,axis1) of the first rotation and the
  // (angle2, axis2) of the second one (if needed)

  int axis1,axis2,angle1,angle2,nrot;
  nrot=2;
  if (new_up == U)
    {
      nrot=1;
      axis1=U;
      switch(new_right)
        {
        case R : angle1=0;  break;
        case B : angle1=1;  break;
        case L : angle1=2;  break;
        case F : angle1=3;  break;
        default : 
         std::ostringstream error_stream;
         error_stream << "new_right is "
                      << new_right << " (" << Direct_string[new_right] 
                      << ") "
                      << "it should be R, B, L, or F" << std::endl;
         throw OomphLibError(error_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
        }
    }
  
  if (new_up == D)
    {
      switch(new_right)
        {
        case R : 
          nrot=1;
          axis1=R; angle1=2;
          break;
        case B :
          axis1=R; angle1=2;
          axis2=U; angle2=1;
          break;
        case L:
          nrot=1;
          axis1=F; angle1=2;
          break;
        case F :
          axis1=R; angle1=2;
          axis2=U; angle2=3;
          break;
        default : 
         std::ostringstream error_stream;
         error_stream << "new_right is "
                      << new_right << " (" << Direct_string[new_right] 
                      << ") "
                      << "it should be R, B, L, or F" << std::endl;
         throw OomphLibError(error_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
         
        }
    }
  
  if (new_up == R)
    {
      switch(new_right)
        {
        case D :
          nrot=1;
          axis1=F; angle1=3;
          break;
        case B :
          axis1=F; angle1=3;
          axis2=R; angle2=1;
          break;
        case U:
          axis1=F; angle1=1;
          axis2=U; angle2=2;
          break;
        case F:
          axis1=F; angle1=3;
          axis2=R; angle2=3;
          break;
        default :
         std::ostringstream error_stream;
         error_stream << "new_right is "
                      << new_right << " (" << Direct_string[new_right] 
                      << ") "
                      << "it should be D, B, U, or F" << std::endl;
         throw OomphLibError(error_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION); 
        }
    }

  if (new_up == L)
    {
      switch(new_right)
        {
        case D :
          axis1=F; angle1=1;
          axis2=R; angle2=2;
          break;
        case B :
          axis1=F; angle1=1;
          axis2=R; angle2=3;
          break;
        case U :
          nrot=1;
          axis1=F; angle1=1;
          break;
        case F :
          axis1=F; angle1=1;
          axis2=R; angle2=1;
          break;
        default :
         std::ostringstream error_stream;
         error_stream << "new_right is "
                      << new_right << " (" << Direct_string[new_right] 
                      << ") "
                      << "it should be D, B, U, or F" << std::endl;
         throw OomphLibError(error_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION); 
        }
    }
  
  if (new_up == F)
      {
        switch(new_right)
          {
          case R :
            nrot=1;
            axis1=R; angle1=1;
            break;
          case L :
            axis1=R; angle1=1;
            axis2=F; angle2=2;
            break;
          case U:
            axis1=R; angle1=1;
            axis2=F; angle2=1;
            break;
          case D:
            axis1=R; angle1=1;
            axis2=F; angle2=3;
            break;
          default : 
           std::ostringstream error_stream;
           error_stream << "new_right is "
                        << new_right << " (" << Direct_string[new_right] 
                        << ") "
                        << "it should be R, L, U, or D" << std::endl;
           throw OomphLibError(error_stream.str(),
                               OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION); 
          }
      }
  if (new_up == B)
    {
      switch(new_right)
        {
        case R :
          nrot=1;
          axis1=R; angle1=3;
          break;
        case L :
          axis1=R; angle1=3;
          axis2=F; angle2=2;
          break;
        case U:
          axis1=R; angle1=3;
          axis2=F; angle2=1;
          break;
        case D:
          axis1=R; angle1=3;
          axis2=F; angle2=3;
          break;
        default : 
         std::ostringstream error_stream;
         error_stream << "new_right is "
                      << new_right << " (" << Direct_string[new_right] 
                      << ") "
                      << "it should be R, L, U, or D" << std::endl;
         throw OomphLibError(error_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION); 
        }
    }
  
  Vector<int> vect_new_dir(3);
  DenseMatrix<int> mat1(3);
  DenseMatrix<int> mat2(3);
  DenseMatrix<int> mat3(3);
  
  // Then we build the first rotation matrix
  construct_rotation_matrix(axis1,angle1,mat1);
  
  // If needed we build the second one
  if (nrot==2)
   {
    // And we make the composition of the two rotations
    // which is stored in Mat3
    construct_rotation_matrix(axis2,angle2,mat2);
    mult_mat_mat(mat2,mat1,mat3);
   }
  else
   {
    // Else we just copy Mat1 into Mat3
    for(int i=0;i<3;i++)
     {
      for(int j=0;j<3;j++)
       {
        mat3(i,j)=mat1(i,j);
       }
     }
   }    
  
  // Rotate
  mult_mat_vect(mat3,dir,vect_new_dir);

  // Return the corresponding vector
  return vect_new_dir;  
}




//====================================================================
/// Setup static data for OcTree 
//====================================================================
void OcTree::setup_static_data()
{
 using namespace OcTreeNames;
 
 
#ifdef PARANOID
 if (Tree::OMEGA!=OcTree::OMEGA)
  {
   std::ostringstream error_stream;
   error_stream 
    << "Inconsistent enumeration!  \n  Tree::OMEGA=" << Tree::OMEGA
                               << "\nOcTree::OMEGA=" << OcTree::OMEGA 
    << std::endl;
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

// Private static data
//--------------------

  // Set flag to indicate that static data has been setup
  Static_data_has_been_setup=true;

 // Initialisation of cosi and sini used in rotation matrix
 Cosi.resize(4);
 Sini.resize(4);
 
 Cosi[0]=1;
 Cosi[1]=0;
 Cosi[2]=-1;
 Cosi[3]=0;

 Sini[0]=0;
 Sini[1]=1;
 Sini[2]=0;
 Sini[3]=-1;


 // Build direction/octant adjacency scheme
 // Is_adjacent(i_face,j_octant): 
 // Is face adjacent to octant?
 // See the table in Samet book
 Is_adjacent.resize(27,27);
 Is_adjacent(L,LDB)=true;
 Is_adjacent(R,LDB)=false;
 Is_adjacent(D,LDB)=true;
 Is_adjacent(U,LDB)=false;
 Is_adjacent(B,LDB)=true;
 Is_adjacent(F,LDB)=false;

 Is_adjacent(L,LDF)=true;
 Is_adjacent(R,LDF)=false;
 Is_adjacent(D,LDF)=true;
 Is_adjacent(U,LDF)=false;
 Is_adjacent(B,LDF)=false;
 Is_adjacent(F,LDF)=true;

 Is_adjacent(L,LUB)=true;
 Is_adjacent(R,LUB)=false;
 Is_adjacent(D,LUB)=false;
 Is_adjacent(U,LUB)=true;
 Is_adjacent(B,LUB)=true;
 Is_adjacent(F,LUB)=false;

 Is_adjacent(L,LUF)=true;
 Is_adjacent(R,LUF)=false;
 Is_adjacent(D,LUF)=false;
 Is_adjacent(U,LUF)=true;
 Is_adjacent(B,LUF)=false;
 Is_adjacent(F,LUF)=true;

 Is_adjacent(L,RDB)=false;
 Is_adjacent(R,RDB)=true;
 Is_adjacent(D,RDB)=true;
 Is_adjacent(U,RDB)=false;
 Is_adjacent(B,RDB)=true;
 Is_adjacent(F,RDB)=false;

 Is_adjacent(L,RDF)=false;
 Is_adjacent(R,RDF)=true;
 Is_adjacent(D,RDF)=true;
 Is_adjacent(U,RDF)=false;
 Is_adjacent(B,RDF)=false;
 Is_adjacent(F,RDF)=true;

 Is_adjacent(L,RUB)=false;
 Is_adjacent(R,RUB)=true;
 Is_adjacent(D,RUB)=false;
 Is_adjacent(U,RUB)=true;
 Is_adjacent(B,RUB)=true;
 Is_adjacent(F,RUB)=false;

 Is_adjacent(L,RUF)=false;
 Is_adjacent(R,RUF)=true;
 Is_adjacent(D,RUF)=false;
 Is_adjacent(U,RUF)=true;
 Is_adjacent(B,RUF)=false;
 Is_adjacent(F,RUF)=true;


 // Build direction/octant adjacency scheme
 // Is_adjacent(i_edge,j_octant): 
 // Is edge adjacent to octant?
 // See the table in Samet book
 Is_adjacent(LD,LDB)=true;
 Is_adjacent(LD,LDF)=true;
 Is_adjacent(LD,LUB)=false;
 Is_adjacent(LD,LUF)=false;
 Is_adjacent(LD,RDB)=false;
 Is_adjacent(LD,RDF)=false;
 Is_adjacent(LD,RUB)=false;
 Is_adjacent(LD,RUF)=false;

 Is_adjacent(LU,LDB)=false;
 Is_adjacent(LU,LDF)=false;
 Is_adjacent(LU,LUB)=true;
 Is_adjacent(LU,LUF)=true;
 Is_adjacent(LU,RDB)=false;
 Is_adjacent(LU,RDF)=false;
 Is_adjacent(LU,RUB)=false;
 Is_adjacent(LU,RUF)=false;


 Is_adjacent(LB,LDB)=true;
 Is_adjacent(LB,LDF)=false;
 Is_adjacent(LB,LUB)=true;
 Is_adjacent(LB,LUF)=false;
 Is_adjacent(LB,RDB)=false;
 Is_adjacent(LB,RDF)=false;
 Is_adjacent(LB,RUB)=false;
 Is_adjacent(LB,RUF)=false;


 Is_adjacent(LF,LDB)=false;
 Is_adjacent(LF,LDF)=true;
 Is_adjacent(LF,LUB)=false;
 Is_adjacent(LF,LUF)=true;
 Is_adjacent(LF,RDB)=false;
 Is_adjacent(LF,RDF)=false;
 Is_adjacent(LF,RUB)=false;
 Is_adjacent(LF,RUF)=false;


 Is_adjacent(RD,LDB)=false;
 Is_adjacent(RD,LDF)=false;
 Is_adjacent(RD,LUB)=false;
 Is_adjacent(RD,LUF)=false;
 Is_adjacent(RD,RDB)=true;
 Is_adjacent(RD,RDF)=true;
 Is_adjacent(RD,RUB)=false;
 Is_adjacent(RD,RUF)=false;


 Is_adjacent(RU,LDB)=false;
 Is_adjacent(RU,LDF)=false;
 Is_adjacent(RU,LUB)=false;
 Is_adjacent(RU,LUF)=false;
 Is_adjacent(RU,RDB)=false;
 Is_adjacent(RU,RDF)=false;
 Is_adjacent(RU,RUB)=true;
 Is_adjacent(RU,RUF)=true;

 Is_adjacent(RB,LDB)=false;
 Is_adjacent(RB,LDF)=false;
 Is_adjacent(RB,LUB)=false;
 Is_adjacent(RB,LUF)=false;
 Is_adjacent(RB,RDB)=true;
 Is_adjacent(RB,RDF)=false;
 Is_adjacent(RB,RUB)=true;
 Is_adjacent(RB,RUF)=false;

 Is_adjacent(RF,LDB)=false;
 Is_adjacent(RF,LDF)=false;
 Is_adjacent(RF,LUB)=false;
 Is_adjacent(RF,LUF)=false;
 Is_adjacent(RF,RDB)=false;
 Is_adjacent(RF,RDF)=true;
 Is_adjacent(RF,RUB)=false;
 Is_adjacent(RF,RUF)=true;


 Is_adjacent(DB,LDB)=true;
 Is_adjacent(DB,LDF)=false;
 Is_adjacent(DB,LUB)=false;
 Is_adjacent(DB,LUF)=false;
 Is_adjacent(DB,RDB)=true;
 Is_adjacent(DB,RDF)=false;
 Is_adjacent(DB,RUB)=false;
 Is_adjacent(DB,RUF)=false;


 Is_adjacent(DF,LDB)=false;
 Is_adjacent(DF,LDF)=true;
 Is_adjacent(DF,LUB)=false;
 Is_adjacent(DF,LUF)=false;
 Is_adjacent(DF,RDB)=false;
 Is_adjacent(DF,RDF)=true;
 Is_adjacent(DF,RUB)=false;
 Is_adjacent(DF,RUF)=false;

 Is_adjacent(UB,LDB)=false;
 Is_adjacent(UB,LDF)=false;
 Is_adjacent(UB,LUB)=true;
 Is_adjacent(UB,LUF)=false;
 Is_adjacent(UB,RDB)=false;
 Is_adjacent(UB,RDF)=false;
 Is_adjacent(UB,RUB)=true;
 Is_adjacent(UB,RUF)=false;

 Is_adjacent(UF,LDB)=false;
 Is_adjacent(UF,LDF)=false;
 Is_adjacent(UF,LUB)=false;
 Is_adjacent(UF,LUF)=true;
 Is_adjacent(UF,RDB)=false;
 Is_adjacent(UF,RDF)=false;
 Is_adjacent(UF,RUB)=false;
 Is_adjacent(UF,RUF)=true;



 // Common face of various octants from Samet's book
 Common_face.resize(27,27);
 Common_face(LDB,LDB)=OMEGA;
 Common_face(LDB,LDF)=OMEGA;
 Common_face(LDB,LUB)=OMEGA;
 Common_face(LDB,LUF)=L;
 Common_face(LDB,RDB)=OMEGA;
 Common_face(LDB,RDF)=D;
 Common_face(LDB,RUB)=B;
 Common_face(LDB,RUF)=OMEGA;

 Common_face(LDF,LDB)=OMEGA;
 Common_face(LDF,LDF)=OMEGA;
 Common_face(LDF,LUB)=L;
 Common_face(LDF,LUF)=OMEGA;
 Common_face(LDF,RDB)=D;
 Common_face(LDF,RDF)=OMEGA;
 Common_face(LDF,RUB)=OMEGA;
 Common_face(LDF,RUF)=F;

 Common_face(LUB,LDB)=OMEGA;
 Common_face(LUB,LDF)=L;
 Common_face(LUB,LUB)=OMEGA;
 Common_face(LUB,LUF)=OMEGA;
 Common_face(LUB,RDB)=B;
 Common_face(LUB,RDF)=OMEGA;
 Common_face(LUB,RUB)=OMEGA;
 Common_face(LUB,RUF)=U;

 Common_face(LUF,LDB)=L;
 Common_face(LUF,LDF)=OMEGA;
 Common_face(LUF,LUB)=OMEGA;
 Common_face(LUF,LUF)=OMEGA;
 Common_face(LUF,RDB)=OMEGA;
 Common_face(LUF,RDF)=F;
 Common_face(LUF,RUB)=U;
 Common_face(LUF,RUF)=OMEGA;

 Common_face(RDB,LDB)=OMEGA;
 Common_face(RDB,LDF)=D;
 Common_face(RDB,LUB)=B;
 Common_face(RDB,LUF)=OMEGA;
 Common_face(RDB,RDB)=OMEGA;
 Common_face(RDB,RDF)=OMEGA;
 Common_face(RDB,RUB)=OMEGA;
 Common_face(RDB,RUF)=R;

 Common_face(RDF,LDB)=D;
 Common_face(RDF,LDF)=OMEGA;
 Common_face(RDF,LUB)=OMEGA;
 Common_face(RDF,LUF)=F;
 Common_face(RDF,RDB)=OMEGA;
 Common_face(RDF,RDF)=OMEGA;
 Common_face(RDF,RUB)=R;
 Common_face(RDF,RUF)=OMEGA;

 Common_face(RUB,LDB)=B;
 Common_face(RUB,LDF)=OMEGA;
 Common_face(RUB,LUB)=OMEGA;
 Common_face(RUB,LUF)=U;
 Common_face(RUB,RDB)=OMEGA;
 Common_face(RUB,RDF)=R;
 Common_face(RUB,RUB)=OMEGA;
 Common_face(RUB,RUF)=OMEGA;

 Common_face(RUF,LDB)=OMEGA;
 Common_face(RUF,LDF)=F;
 Common_face(RUF,LUB)=U;
 Common_face(RUF,LUF)=OMEGA;
 Common_face(RUF,RDB)=R;
 Common_face(RUF,RDF)=OMEGA;
 Common_face(RUF,RUB)=OMEGA;
 Common_face(RUF,RUF)=OMEGA;



 // Common face of various edges/octants from Samet's book
 Common_face(LD,LDB)=OMEGA;
 Common_face(LD,LDF)=OMEGA;
 Common_face(LD,LUB)=L;
 Common_face(LD,LUF)=L;
 Common_face(LD,RDB)=D;
 Common_face(LD,RDF)=D;
 Common_face(LD,RUB)=OMEGA;
 Common_face(LD,RUF)=OMEGA;

 Common_face(LU,LDB)=L;
 Common_face(LU,LDF)=L;
 Common_face(LU,LUB)=OMEGA;
 Common_face(LU,LUF)=OMEGA;
 Common_face(LU,RDB)=OMEGA;
 Common_face(LU,RDF)=OMEGA;
 Common_face(LU,RUB)=U;
 Common_face(LU,RUF)=U;

 Common_face(LB,LDB)=OMEGA;
 Common_face(LB,LDF)=L;
 Common_face(LB,LUB)=OMEGA;
 Common_face(LB,LUF)=L;
 Common_face(LB,RDB)=B;
 Common_face(LB,RDF)=OMEGA;
 Common_face(LB,RUB)=B;
 Common_face(LB,RUF)=OMEGA;

 Common_face(LF,LDB)=L;
 Common_face(LF,LDF)=OMEGA;
 Common_face(LF,LUB)=L;
 Common_face(LF,LUF)=OMEGA;
 Common_face(LF,RDB)=OMEGA;
 Common_face(LF,RDF)=F;
 Common_face(LF,RUB)=OMEGA;
 Common_face(LF,RUF)=F;

 Common_face(RD,LDB)=D;
 Common_face(RD,LDF)=D;
 Common_face(RD,LUB)=OMEGA;
 Common_face(RD,LUF)=OMEGA;
 Common_face(RD,RDB)=OMEGA;
 Common_face(RD,RDF)=OMEGA;
 Common_face(RD,RUB)=R;
 Common_face(RD,RUF)=R;

 Common_face(RU,LDB)=OMEGA;
 Common_face(RU,LDF)=OMEGA;
 Common_face(RU,LUB)=U;
 Common_face(RU,LUF)=U;
 Common_face(RU,RDB)=R;
 Common_face(RU,RDF)=R;
 Common_face(RU,RUB)=OMEGA;
 Common_face(RU,RUF)=OMEGA;

 Common_face(RB,LDB)=B;
 Common_face(RB,LDF)=OMEGA;
 Common_face(RB,LUB)=B;
 Common_face(RB,LUF)=OMEGA;
 Common_face(RB,RDB)=OMEGA;
 Common_face(RB,RDF)=R;
 Common_face(RB,RUB)=OMEGA;
 Common_face(RB,RUF)=R;

 Common_face(RF,LDB)=OMEGA;
 Common_face(RF,LDF)=F;
 Common_face(RF,LUB)=OMEGA;
 Common_face(RF,LUF)=F;
 Common_face(RF,RDB)=R;
 Common_face(RF,RDF)=OMEGA;
 Common_face(RF,RUB)=R;
 Common_face(RF,RUF)=OMEGA;


 Common_face(DB,LDB)=OMEGA;
 Common_face(DB,LDF)=D;
 Common_face(DB,LUB)=B;
 Common_face(DB,LUF)=OMEGA;
 Common_face(DB,RDB)=OMEGA;
 Common_face(DB,RDF)=D;
 Common_face(DB,RUB)=B;
 Common_face(DB,RUF)=OMEGA;

 Common_face(DF,LDB)=D;
 Common_face(DF,LDF)=OMEGA;
 Common_face(DF,LUB)=OMEGA;
 Common_face(DF,LUF)=F;
 Common_face(DF,RDB)=D;
 Common_face(DF,RDF)=OMEGA;
 Common_face(DF,RUB)=OMEGA;
 Common_face(DF,RUF)=F;

 Common_face(UB,LDB)=B;
 Common_face(UB,LDF)=OMEGA;
 Common_face(UB,LUB)=OMEGA;
 Common_face(UB,LUF)=U;
 Common_face(UB,RDB)=B;
 Common_face(UB,RDF)=OMEGA;
 Common_face(UB,RUB)=OMEGA;
 Common_face(UB,RUF)=U;

 Common_face(UF,LDB)=OMEGA;
 Common_face(UF,LDF)=F;
 Common_face(UF,LUB)=U;
 Common_face(UF,LUF)=OMEGA;
 Common_face(UF,RDB)=OMEGA;
 Common_face(UF,RDF)=F;
 Common_face(UF,RUB)=U;
 Common_face(UF,RUF)=OMEGA;


 // Tecplot colours for neighbours in various directions
 Colour.resize(27);
 Colour[R]="CYAN";
 Colour[L]="RED";
 Colour[U]="GREEN";
 Colour[D]="BLUE";
 Colour[F]="CUSTOM3";
 Colour[B]="PURPLE";
 Colour[OMEGA]="YELLOW";

 Colour[LB]="BLUE";
 Colour[RB]="BLUE";
 Colour[DB]="BLUE";
 Colour[UB]="BLUE";

 Colour[LD]="GREEN";
 Colour[RD]="GREEN";
 Colour[LU]="GREEN";
 Colour[RU]="GREEN";


 Colour[LF]="RED";
 Colour[RF]="RED";
 Colour[DF]="RED";
 Colour[UF]="RED";

 Colour[OMEGA]="YELLOW";



 // Reflection scheme:
 // Reflect(direction,octant): Get mirror of octant in direction 
 // Table in Samet book as well
 Reflect.resize(27,27);

 Reflect(L,LDB)=RDB;
 Reflect(R,LDB)=RDB;
 Reflect(D,LDB)=LUB;
 Reflect(U,LDB)=LUB;
 Reflect(B,LDB)=LDF;
 Reflect(F,LDB)=LDF;

 Reflect(L,LDF)=RDF;
 Reflect(R,LDF)=RDF;
 Reflect(D,LDF)=LUF;
 Reflect(U,LDF)=LUF;
 Reflect(B,LDF)=LDB;
 Reflect(F,LDF)=LDB;

 Reflect(L,LUB)=RUB;
 Reflect(R,LUB)=RUB;
 Reflect(D,LUB)=LDB;
 Reflect(U,LUB)=LDB;
 Reflect(B,LUB)=LUF;
 Reflect(F,LUB)=LUF;
 
 Reflect(L,LUF)=RUF;
 Reflect(R,LUF)=RUF;
 Reflect(D,LUF)=LDF;
 Reflect(U,LUF)=LDF;
 Reflect(B,LUF)=LUB;
 Reflect(F,LUF)=LUB;

 Reflect(L,RDB)=LDB;
 Reflect(R,RDB)=LDB;
 Reflect(D,RDB)=RUB;
 Reflect(U,RDB)=RUB;
 Reflect(B,RDB)=RDF;
 Reflect(F,RDB)=RDF;

 Reflect(L,RDF)=LDF;
 Reflect(R,RDF)=LDF;
 Reflect(D,RDF)=RUF;
 Reflect(U,RDF)=RUF;
 Reflect(B,RDF)=RDB;
 Reflect(F,RDF)=RDB;

 Reflect(L,RUB)=LUB;
 Reflect(R,RUB)=LUB;
 Reflect(D,RUB)=RDB;
 Reflect(U,RUB)=RDB;
 Reflect(B,RUB)=RUF;
 Reflect(F,RUB)=RUF;
 
 Reflect(L,RUF)=LUF;
 Reflect(R,RUF)=LUF;
 Reflect(D,RUF)=RDF;
 Reflect(U,RUF)=RDF;
 Reflect(B,RUF)=RUB;
 Reflect(F,RUF)=RUB;


 // Reflection scheme:
 // Reflect(direction (edge) ,octant): Get mirror of edge in direction 
 // Table in Samet book as well

 Reflect(LD,LDB)=RUB;
 Reflect(LU,LDB)=RUB;
 Reflect(RD,LDB)=RUB;
 Reflect(RU,LDB)=RUB;

 Reflect(LD,LDF)=RUF;
 Reflect(LU,LDF)=RUF;
 Reflect(RD,LDF)=RUF;
 Reflect(RU,LDF)=RUF;

 Reflect(LD,LUB)=RDB;
 Reflect(LU,LUB)=RDB;
 Reflect(RD,LUB)=RDB;
 Reflect(RU,LUB)=RDB;

 Reflect(LD,LUF)=RDF;
 Reflect(LU,LUF)=RDF;
 Reflect(RD,LUF)=RDF;
 Reflect(RU,LUF)=RDF;


 Reflect(LD,RDB)=LUB;
 Reflect(LU,RDB)=LUB;
 Reflect(RD,RDB)=LUB;
 Reflect(RU,RDB)=LUB;

 Reflect(LD,RDF)=LUF;
 Reflect(LU,RDF)=LUF;
 Reflect(RD,RDF)=LUF;
 Reflect(RU,RDF)=LUF;

 Reflect(LD,RUB)=LDB;
 Reflect(LU,RUB)=LDB;
 Reflect(RD,RUB)=LDB;
 Reflect(RU,RUB)=LDB;

 Reflect(LD,RUF)=LDF;
 Reflect(LU,RUF)=LDF;
 Reflect(RD,RUF)=LDF;
 Reflect(RU,RUF)=LDF;



 Reflect(LB,LDB)=RDF;
 Reflect(LF,LDB)=RDF;
 Reflect(RB,LDB)=RDF;
 Reflect(RF,LDB)=RDF;

 Reflect(LB,LDF)=RDB;
 Reflect(LF,LDF)=RDB;
 Reflect(RB,LDF)=RDB;
 Reflect(RF,LDF)=RDB;

 Reflect(LB,LUB)=RUF;
 Reflect(LF,LUB)=RUF;
 Reflect(RB,LUB)=RUF;
 Reflect(RF,LUB)=RUF;

 Reflect(LB,LUF)=RUB;
 Reflect(LF,LUF)=RUB;
 Reflect(RB,LUF)=RUB;
 Reflect(RF,LUF)=RUB;

 Reflect(LB,RDB)=LDF;
 Reflect(LF,RDB)=LDF;
 Reflect(RB,RDB)=LDF;
 Reflect(RF,RDB)=LDF;

 Reflect(LB,RDF)=LDB;
 Reflect(LF,RDF)=LDB;
 Reflect(RB,RDF)=LDB;
 Reflect(RF,RDF)=LDB;

 Reflect(LB,RUB)=LUF;
 Reflect(LF,RUB)=LUF;
 Reflect(RB,RUB)=LUF;
 Reflect(RF,RUB)=LUF;

 Reflect(LB,RUF)=LUB;
 Reflect(LF,RUF)=LUB;
 Reflect(RB,RUF)=LUB;
 Reflect(RF,RUF)=LUB;


 Reflect(DB,LDB)=LUF;
 Reflect(DF,LDB)=LUF;
 Reflect(UB,LDB)=LUF;
 Reflect(UF,LDB)=LUF;

 Reflect(DB,LDF)=LUB;
 Reflect(DF,LDF)=LUB;
 Reflect(UB,LDF)=LUB;
 Reflect(UF,LDF)=LUB;

 Reflect(DB,LUB)=LDF;
 Reflect(DF,LUB)=LDF;
 Reflect(UB,LUB)=LDF;
 Reflect(UF,LUB)=LDF;

 Reflect(DB,LUF)=LDB;
 Reflect(DF,LUF)=LDB;
 Reflect(UB,LUF)=LDB;
 Reflect(UF,LUF)=LDB;

 Reflect(DB,RDB)=RUF;
 Reflect(DF,RDB)=RUF;
 Reflect(UB,RDB)=RUF;
 Reflect(UF,RDB)=RUF;

 Reflect(DB,RDF)=RUB;
 Reflect(DF,RDF)=RUB;
 Reflect(UB,RDF)=RUB;
 Reflect(UF,RDF)=RUB;

 Reflect(DB,RUB)=RDF;
 Reflect(DF,RUB)=RDF;
 Reflect(UB,RUB)=RDF;
 Reflect(UF,RUB)=RDF;

 Reflect(DB,RUF)=RDB;
 Reflect(DF,RUF)=RDB;
 Reflect(UB,RUF)=RDB;
 Reflect(UF,RUF)=RDB;



 // S_base(i,direction) etc.:  Initial value/increment for coordinate s[i] on 
 // the face indicated by direction (L/R/U/D/F/B) 
 S_base.resize(3,27);
 S_steplo.resize(3,27);
 S_stephi.resize(3,27);

 S_base(0,L)=-1.0; 
 S_base(1,L)=-1.0; 
 S_base(2,L)=-1.0; 
 S_steplo(0,L)=0.0;
 S_steplo(1,L)=2.0;
 S_steplo(2,L)=0.0;
 S_stephi(0,L)=0.0;
 S_stephi(1,L)=0.0;
 S_stephi(2,L)=2.0;

 S_base(0,R)=1.0;
 S_base(1,R)=-1.0;
 S_base(2,R)=-1.0;
 S_steplo(0,R)=0.0;
 S_steplo(1,R)=2.0;
 S_steplo(2,R)=0.0;
 S_stephi(0,R)=0.0;
 S_stephi(1,R)=0.0;
 S_stephi(2,R)=2.0;

 S_base(0,U)=-1.0;
 S_base(1,U)=1.0;
 S_base(2,U)=-1.0;
 S_steplo(0,U)=2.0;
 S_steplo(1,U)=0.0;
 S_steplo(2,U)=0.0;
 S_stephi(0,U)=0.0;
 S_stephi(1,U)=0.0;
 S_stephi(2,U)=2.0;

 S_base(0,D)=-1.0;
 S_base(1,D)=-1.0;
 S_base(2,D)=-1.0;
 S_steplo(0,D)=2.0;
 S_steplo(1,D)=0.0;
 S_steplo(2,D)=0.0;
 S_stephi(0,D)=0.0;
 S_stephi(1,D)=0.0;
 S_stephi(2,D)=2.0;

 S_base(0,F)=-1.0;
 S_base(1,F)=-1.0;
 S_base(2,F)=1.0;
 S_steplo(0,F)=2.0;
 S_steplo(1,F)=0.0;
 S_steplo(2,F)=0.0;
 S_stephi(0,F)=0.0;
 S_stephi(1,F)=2.0;
 S_stephi(2,F)=0.0;

 S_base(0,B)=-1.0;
 S_base(1,B)=-1.0;
 S_base(2,B)=-1.0;
 S_steplo(0,B)=2.0;
 S_steplo(1,B)=0.0;
 S_steplo(2,B)=0.0;
 S_stephi(0,B)=0.0;
 S_stephi(1,B)=2.0;
 S_stephi(2,B)=0.0;

 
 // Relative to the left/down/back vertex in any (father) octree, the 
 // corresponding vertex in the son specified by \c son_octant has an offset.
 // If we project the son_octant's left/down/back vertex onto the
 // father's face \c face, it is located at the in-face coordinate 
 //  \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
 //  \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
 S_directlo.resize(27,27);
 S_directhi.resize(27,27);

 S_directlo(L,LDB)=0.0;
 S_directlo(R,LDB)=0.0;
 S_directlo(U,LDB)=0.0;
 S_directlo(D,LDB)=0.0;
 S_directlo(F,LDB)=0.0;
 S_directlo(B,LDB)=0.0;

 S_directlo(L,LDF)=0.0;
 S_directlo(R,LDF)=0.0;
 S_directlo(U,LDF)=0.0;
 S_directlo(D,LDF)=0.0;
 S_directlo(F,LDF)=0.0;
 S_directlo(B,LDF)=0.0;

 S_directlo(L,LUB)=1.0;
 S_directlo(R,LUB)=1.0;
 S_directlo(U,LUB)=0.0;
 S_directlo(D,LUB)=0.0;
 S_directlo(F,LUB)=0.0;
 S_directlo(B,LUB)=0.0;

 S_directlo(L,LUF)=1.0;
 S_directlo(R,LUF)=1.0;
 S_directlo(U,LUF)=0.0;
 S_directlo(D,LUF)=0.0;
 S_directlo(F,LUF)=0.0;
 S_directlo(B,LUF)=0.0;

 S_directlo(L,RDB)=0.0;
 S_directlo(R,RDB)=0.0;
 S_directlo(U,RDB)=1.0;
 S_directlo(D,RDB)=1.0;
 S_directlo(F,RDB)=1.0;
 S_directlo(B,RDB)=1.0;

 S_directlo(L,RDF)=0.0;
 S_directlo(R,RDF)=0.0;
 S_directlo(U,RDF)=1.0;
 S_directlo(D,RDF)=1.0;
 S_directlo(F,RDF)=1.0;
 S_directlo(B,RDF)=1.0;

 S_directlo(L,RUB)=1.0;
 S_directlo(R,RUB)=1.0;
 S_directlo(U,RUB)=1.0;
 S_directlo(D,RUB)=1.0;
 S_directlo(F,RUB)=1.0;
 S_directlo(B,RUB)=1.0;

 S_directlo(L,RUF)=1.0;
 S_directlo(R,RUF)=1.0;
 S_directlo(U,RUF)=1.0;
 S_directlo(D,RUF)=1.0;
 S_directlo(F,RUF)=1.0;
 S_directlo(B,RUF)=1.0;




 S_directhi(L,LDB)=0.0;
 S_directhi(R,LDB)=0.0;
 S_directhi(U,LDB)=0.0;
 S_directhi(D,LDB)=0.0;
 S_directhi(F,LDB)=0.0;
 S_directhi(B,LDB)=0.0;

 S_directhi(L,LDF)=1.0;
 S_directhi(R,LDF)=1.0;
 S_directhi(U,LDF)=1.0;
 S_directhi(D,LDF)=1.0;
 S_directhi(F,LDF)=0.0;
 S_directhi(B,LDF)=0.0;

 S_directhi(L,LUB)=0.0;
 S_directhi(R,LUB)=0.0;
 S_directhi(U,LUB)=0.0;
 S_directhi(D,LUB)=0.0;
 S_directhi(F,LUB)=1.0;
 S_directhi(B,LUB)=1.0;

 S_directhi(L,LUF)=1.0;
 S_directhi(R,LUF)=1.0;
 S_directhi(U,LUF)=1.0;
 S_directhi(D,LUF)=1.0;
 S_directhi(F,LUF)=1.0;
 S_directhi(B,LUF)=1.0;

 S_directhi(L,RDB)=0.0;
 S_directhi(R,RDB)=0.0;
 S_directhi(U,RDB)=0.0;
 S_directhi(D,RDB)=0.0;
 S_directhi(F,RDB)=0.0;
 S_directhi(B,RDB)=0.0;

 S_directhi(L,RDF)=1.0;
 S_directhi(R,RDF)=1.0;
 S_directhi(U,RDF)=1.0;
 S_directhi(D,RDF)=1.0;
 S_directhi(F,RDF)=0.0;
 S_directhi(B,RDF)=0.0;

 S_directhi(L,RUB)=0.0;
 S_directhi(R,RUB)=0.0;
 S_directhi(U,RUB)=0.0;
 S_directhi(D,RUB)=0.0;
 S_directhi(F,RUB)=1.0;
 S_directhi(B,RUB)=1.0;

 S_directhi(L,RUF)=1.0;
 S_directhi(R,RUF)=1.0;
 S_directhi(U,RUF)=1.0;
 S_directhi(D,RUF)=1.0;
 S_directhi(F,RUF)=1.0;
 S_directhi(B,RUF)=1.0;







 // S_base_edge(i,direction):  Initial value/increment for coordinate s[i] on 
 // the edge indicated by direction (LB,RB,...)
 S_base_edge.resize(3,27);
 S_step_edge.resize(3,27);

 S_base_edge(0,LB)=-1.0; 
 S_base_edge(1,LB)=-1.0; 
 S_base_edge(2,LB)=-1.0; 
 S_step_edge(0,LB)=0.0;
 S_step_edge(1,LB)=2.0;
 S_step_edge(2,LB)=0.0;

 S_base_edge(0,RB)= 1.0; 
 S_base_edge(1,RB)=-1.0; 
 S_base_edge(2,RB)=-1.0; 
 S_step_edge(0,RB)=0.0;
 S_step_edge(1,RB)=2.0;
 S_step_edge(2,RB)=0.0;

 S_base_edge(0,DB)=-1.0; 
 S_base_edge(1,DB)=-1.0; 
 S_base_edge(2,DB)=-1.0; 
 S_step_edge(0,DB)=2.0;
 S_step_edge(1,DB)=0.0;
 S_step_edge(2,DB)=0.0;

 S_base_edge(0,UB)=-1.0; 
 S_base_edge(1,UB)= 1.0; 
 S_base_edge(2,UB)=-1.0; 
 S_step_edge(0,UB)=2.0;
 S_step_edge(1,UB)=0.0;
 S_step_edge(2,UB)=0.0;

 S_base_edge(0,LD)=-1.0; 
 S_base_edge(1,LD)=-1.0; 
 S_base_edge(2,LD)=-1.0; 
 S_step_edge(0,LD)=0.0;
 S_step_edge(1,LD)=0.0;
 S_step_edge(2,LD)=2.0;

 S_base_edge(0,RD)= 1.0; 
 S_base_edge(1,RD)=-1.0; 
 S_base_edge(2,RD)=-1.0; 
 S_step_edge(0,RD)=0.0;
 S_step_edge(1,RD)=0.0;
 S_step_edge(2,RD)=2.0;

 S_base_edge(0,LU)=-1.0; 
 S_base_edge(1,LU)= 1.0; 
 S_base_edge(2,LU)=-1.0; 
 S_step_edge(0,LU)=0.0;
 S_step_edge(1,LU)=0.0;
 S_step_edge(2,LU)=2.0;

 S_base_edge(0,RU)= 1.0; 
 S_base_edge(1,RU)= 1.0; 
 S_base_edge(2,RU)=-1.0; 
 S_step_edge(0,RU)=0.0;
 S_step_edge(1,RU)=0.0;
 S_step_edge(2,RU)=2.0;



 S_base_edge(0,LF)=-1.0; 
 S_base_edge(1,LF)=-1.0; 
 S_base_edge(2,LF)= 1.0; 
 S_step_edge(0,LF)=0.0;
 S_step_edge(1,LF)=2.0;
 S_step_edge(2,LF)=0.0;

 S_base_edge(0,RF)= 1.0; 
 S_base_edge(1,RF)=-1.0; 
 S_base_edge(2,RF)= 1.0; 
 S_step_edge(0,RF)=0.0;
 S_step_edge(1,RF)=2.0;
 S_step_edge(2,RF)=0.0;

 S_base_edge(0,DF)=-1.0; 
 S_base_edge(1,DF)=-1.0; 
 S_base_edge(2,DF)= 1.0; 
 S_step_edge(0,DF)=2.0;
 S_step_edge(1,DF)=0.0;
 S_step_edge(2,DF)=0.0;

 S_base_edge(0,UF)=-1.0; 
 S_base_edge(1,UF)= 1.0; 
 S_base_edge(2,UF)= 1.0; 
 S_step_edge(0,UF)=2.0;
 S_step_edge(1,UF)=0.0;
 S_step_edge(2,UF)=0.0;




 // Relative to the left/down/back vertex in any (father) octree, the 
 // corresponding vertex in the son specified by \c son_octant has an offset.
 // If we project the son_octant's left/down/back vertex onto the
 // father's edge \c edge, it is located at the in-face coordinate 
 //  \c s_lo = h/2 \c S_direct_edge(edge,son_octant). 
 S_direct_edge.resize(27,27);
 S_direct_edge(LB,LDB)=0.0;
 S_direct_edge(RB,LDB)=0.0;
 S_direct_edge(DB,LDB)=0.0;
 S_direct_edge(UB,LDB)=0.0;
 S_direct_edge(LD,LDB)=0.0;
 S_direct_edge(RD,LDB)=0.0;
 S_direct_edge(LU,LDB)=0.0;
 S_direct_edge(RU,LDB)=0.0;
 S_direct_edge(LF,LDB)=0.0;
 S_direct_edge(RF,LDB)=0.0;
 S_direct_edge(DF,LDB)=0.0;
 S_direct_edge(UF,LDB)=0.0;

 S_direct_edge(LB,RDB)=0.0;
 S_direct_edge(RB,RDB)=0.0;
 S_direct_edge(DB,RDB)=1.0;
 S_direct_edge(UB,RDB)=1.0;
 S_direct_edge(LD,RDB)=0.0;
 S_direct_edge(RD,RDB)=0.0;
 S_direct_edge(LU,RDB)=0.0;
 S_direct_edge(RU,RDB)=0.0;
 S_direct_edge(LF,RDB)=0.0;
 S_direct_edge(RF,RDB)=0.0;
 S_direct_edge(DF,RDB)=1.0;
 S_direct_edge(UF,RDB)=1.0;

 S_direct_edge(LB,LUB)=1.0;
 S_direct_edge(RB,LUB)=1.0;
 S_direct_edge(DB,LUB)=0.0;
 S_direct_edge(UB,LUB)=0.0;
 S_direct_edge(LD,LUB)=0.0;
 S_direct_edge(RD,LUB)=0.0;
 S_direct_edge(LU,LUB)=0.0;
 S_direct_edge(RU,LUB)=0.0;
 S_direct_edge(LF,LUB)=1.0;
 S_direct_edge(RF,LUB)=1.0;
 S_direct_edge(DF,LUB)=0.0;
 S_direct_edge(UF,LUB)=0.0;

 S_direct_edge(LB,RUB)=1.0;
 S_direct_edge(RB,RUB)=1.0;
 S_direct_edge(DB,RUB)=1.0;
 S_direct_edge(UB,RUB)=1.0;
 S_direct_edge(LD,RUB)=0.0;
 S_direct_edge(RD,RUB)=0.0;
 S_direct_edge(LU,RUB)=0.0;
 S_direct_edge(RU,RUB)=0.0;
 S_direct_edge(LF,RUB)=1.0;
 S_direct_edge(RF,RUB)=1.0;
 S_direct_edge(DF,RUB)=1.0;
 S_direct_edge(UF,RUB)=1.0;


 S_direct_edge(LB,LDF)=0.0;
 S_direct_edge(RB,LDF)=0.0;
 S_direct_edge(DB,LDF)=0.0;
 S_direct_edge(UB,LDF)=0.0;
 S_direct_edge(LD,LDF)=1.0;
 S_direct_edge(RD,LDF)=1.0;
 S_direct_edge(LU,LDF)=1.0;
 S_direct_edge(RU,LDF)=1.0;
 S_direct_edge(LF,LDF)=0.0;
 S_direct_edge(RF,LDF)=0.0;
 S_direct_edge(DF,LDF)=0.0;
 S_direct_edge(UF,LDF)=0.0;

 S_direct_edge(LB,RDF)=0.0;
 S_direct_edge(RB,RDF)=0.0;
 S_direct_edge(DB,RDF)=1.0;
 S_direct_edge(UB,RDF)=1.0;
 S_direct_edge(LD,RDF)=1.0;
 S_direct_edge(RD,RDF)=1.0;
 S_direct_edge(LU,RDF)=1.0;
 S_direct_edge(RU,RDF)=1.0;
 S_direct_edge(LF,RDF)=0.0;
 S_direct_edge(RF,RDF)=0.0;
 S_direct_edge(DF,RDF)=1.0;
 S_direct_edge(UF,RDF)=1.0;

 S_direct_edge(LB,LUF)=1.0;
 S_direct_edge(RB,LUF)=1.0;
 S_direct_edge(DB,LUF)=0.0;
 S_direct_edge(UB,LUF)=0.0;
 S_direct_edge(LD,LUF)=1.0;
 S_direct_edge(RD,LUF)=1.0;
 S_direct_edge(LU,LUF)=1.0;
 S_direct_edge(RU,LUF)=1.0;
 S_direct_edge(LF,LUF)=1.0;
 S_direct_edge(RF,LUF)=1.0;
 S_direct_edge(DF,LUF)=0.0;
 S_direct_edge(UF,LUF)=0.0;

 S_direct_edge(LB,RUF)=1.0;
 S_direct_edge(RB,RUF)=1.0;
 S_direct_edge(DB,RUF)=1.0;
 S_direct_edge(UB,RUF)=1.0;
 S_direct_edge(LD,RUF)=1.0;
 S_direct_edge(RD,RUF)=1.0;
 S_direct_edge(LU,RUF)=1.0;
 S_direct_edge(RU,RUF)=1.0;
 S_direct_edge(LF,RUF)=1.0;
 S_direct_edge(RF,RUF)=1.0;
 S_direct_edge(DF,RUF)=1.0;
 S_direct_edge(UF,RUF)=1.0;



//Public static data:
//-------------------

 // Translate (enumerated) directions into strings
 Direct_string.resize(27);
 Direct_string[LDB]="LDB";
 Direct_string[LDF]="LDF";
 Direct_string[LUB]="LUB";
 Direct_string[LUF]="LUF";
 Direct_string[RDB]="RDB";
 Direct_string[RDF]="RDF";
 Direct_string[RUB]="RUB";
 Direct_string[RUF]="RUF";


 Direct_string[L]="L";
 Direct_string[R]="R";
 Direct_string[U]="U";
 Direct_string[D]="D";
 Direct_string[F]="F";
 Direct_string[B]="B";
 
 Direct_string[LU]="LU";
 Direct_string[LD]="LD";
 Direct_string[LF]="LF";
 Direct_string[LB]="LB";
 Direct_string[RU]="RU";
 Direct_string[RD]="RD";
 Direct_string[RF]="RF";
 Direct_string[RB]="RB";
 Direct_string[UF]="UF";
 Direct_string[UB]="UB";
 Direct_string[DF]="DF";
 Direct_string[DB]="DB";

 Direct_string[OMEGA]="OMEGA";

 
 // Get opposite face, e.g. Reflect_face(L)=R 
 Reflect_face.resize(27);
 Reflect_face[L]=R;
 Reflect_face[R]=L;
 Reflect_face[U]=D;
 Reflect_face[D]=U;
 Reflect_face[B]=F;
 Reflect_face[F]=B;

 // Get opposite edge, e.g. Reflect_edge(DB)=UF 
 Reflect_edge.resize(27);
 Reflect_edge[LB]=RF;
 Reflect_edge[RB]=LF;
 Reflect_edge[DB]=UF;
 Reflect_edge[UB]=DF;

 Reflect_edge[LD]=RU;
 Reflect_edge[RD]=LU;
 Reflect_edge[LU]=RD;
 Reflect_edge[RU]=LD;

 Reflect_edge[LF]=RB;
 Reflect_edge[RF]=LB;
 Reflect_edge[DF]=UB;
 Reflect_edge[UF]=DB;

 // Get opposite vertex, e.g. Reflect_vertex(LDB)=RUF
 Reflect_vertex.resize(27);
 Reflect_vertex[LDB]=RUF;
 Reflect_vertex[RUF]=LDB;
 Reflect_vertex[RDB]=LUF;
 Reflect_vertex[LUF]=RDB;
 Reflect_vertex[LUB]=RDF;
 Reflect_vertex[RDF]=LUB;
 Reflect_vertex[RUB]=LDF;
 Reflect_vertex[LDF]=RUB;



 // Vertices at ends of edges
 Vertex_at_end_of_edge.resize(27);

 Vertex_at_end_of_edge[DB].resize(2);
 Vertex_at_end_of_edge[DB][0]=LDB; // Pattern: both other indices
 Vertex_at_end_of_edge[DB][1]=RDB;

 Vertex_at_end_of_edge[UB].resize(2);
 Vertex_at_end_of_edge[UB][0]=LUB; // Pattern: both other indices
 Vertex_at_end_of_edge[UB][1]=RUB;


 Vertex_at_end_of_edge[LB].resize(2);
 Vertex_at_end_of_edge[LB][0]=LUB; // Pattern: both other indices
 Vertex_at_end_of_edge[LB][1]=LDB;

 Vertex_at_end_of_edge[RB].resize(2);
 Vertex_at_end_of_edge[RB][0]=RUB; // Pattern: both other indices
 Vertex_at_end_of_edge[RB][1]=RDB;



 Vertex_at_end_of_edge[LD].resize(2);
 Vertex_at_end_of_edge[LD][0]=LDF; // Pattern: both other indices
 Vertex_at_end_of_edge[LD][1]=LDB;

 Vertex_at_end_of_edge[RD].resize(2);
 Vertex_at_end_of_edge[RD][0]=RDF; // Pattern: both other indices
 Vertex_at_end_of_edge[RD][1]=RDB;

 Vertex_at_end_of_edge[LU].resize(2);
 Vertex_at_end_of_edge[LU][0]=LUF; // Pattern: both other indices
 Vertex_at_end_of_edge[LU][1]=LUB;

 Vertex_at_end_of_edge[RU].resize(2);
 Vertex_at_end_of_edge[RU][0]=RUF; // Pattern: both other indices
 Vertex_at_end_of_edge[RU][1]=RUB;



 Vertex_at_end_of_edge[DF].resize(2);
 Vertex_at_end_of_edge[DF][0]=LDF; // Pattern: both other indices
 Vertex_at_end_of_edge[DF][1]=RDF;

 Vertex_at_end_of_edge[UF].resize(2);
 Vertex_at_end_of_edge[UF][0]=LUF; // Pattern: both other indices
 Vertex_at_end_of_edge[UF][1]=RUF;

 Vertex_at_end_of_edge[LF].resize(2);
 Vertex_at_end_of_edge[LF][0]=LUF; // Pattern: both other indices
 Vertex_at_end_of_edge[LF][1]=LDF;

 Vertex_at_end_of_edge[RF].resize(2);
 Vertex_at_end_of_edge[RF][0]=RUF; // Pattern: both other indices
 Vertex_at_end_of_edge[RF][1]=RDF;




 // Initialisation of the values of Vector_to_direction
 Vector<int> vect(3);
 int elem;

 for(int i=-1; i<2; i++)
   {
     for(int j=-1; j<2; j++)
       {
         for(int k=-1; k<2; k++)
           {
             vect[0]=i; vect[1]=j; vect[2]=k;
             int num_elem=0;
             
             // To put a number on the vector (i,j,k), we assume that that
             // the vector (i+1,j+1,k+1) represents the decomposition
             // of the number of the corresponding direction in base 3.
             num_elem=(i+1)*9+(j+1)*3+(k+1);
          
             //for each number we have the corresponding element
             switch(num_elem)
               {
               case 6  :elem=LUB; break;
               case 24 :elem=RUB; break;
               case 26 :elem=RUF; break;
               case 8  :elem=LUF; break;
               case 0  :elem=LDB; break;
               case 18 :elem=RDB; break;
               case 20 :elem=RDF; break;
               case 2  :elem=LDF; break;
               case 25 :elem=RU; break;
               case 23 :elem=RF; break;
               case 19 :elem=RD; break;
               case 21 :elem=RB; break;
               case 7  :elem=LU; break;
               case 5  :elem=LF; break;
               case 1  :elem=LD; break;
               case 3  :elem=LB; break;
               case 17 :elem=UF; break;
               case 15 :elem=UB; break;
               case 11 :elem=DF; break;
               case 9  :elem=DB; break;
               case 16 :elem=U; break;
               case 10 :elem=D; break;
               case 22 :elem=R; break;
               case 4  :elem=L; break;
               case 14 :elem=F; break;
               case 12 :elem=B; break;
               case 13 :elem=OMEGA; break;
               default : 
                elem=OMEGA; 
                oomph_info 
                 << "there might be a problem with Vector_to_direction" 
                 <<std::endl;
                 break;
               }
             Vector_to_direction[vect]=elem;
           }
       }
   }



 // Initialisation of Direction_to_vector:
 // Translate Octant, face, edge into direction vector using the 
 // value of Direct_string; Direction_to_vector[U]={0,1,0};
 Direction_to_vector.resize(27);
 for(int i=LDB;i<=F;i++)
   {
     Direction_to_vector[i].resize(3);
     // Initialisation to 0;    
     for(int j=0;j<3;j++)
       {
         Direction_to_vector[i][j]=0;
       }
     
     // Use +1 or -1 to indicate the relevant components of
     // the vector Direction 
     for (unsigned j = 0; j < 3; j++)
       {
         if (Direct_string[i].length() > j) 
           {
             switch (Direct_string[i][j])
               {
               case 'R':
                 Direction_to_vector[i][0] = 1;
                 break;
               case 'L':
                 Direction_to_vector[i][0] = -1;
                 break;
               case 'U':
                 Direction_to_vector[i][1] = 1;
                 break;
               case 'D':
                 Direction_to_vector[i][1] = -1;
                 break;
               case 'F':
                 Direction_to_vector[i][2] = 1;
                 break;
               case 'B':
                 Direction_to_vector[i][2] = -1;
                 break;
               default:
                 oomph_info << "Direction Error !!" << std::endl;
               }
           }
       }
   }


 // Setup map that works out required rotations based on
 //-----------------------------------------------------
 // adjacent edge vertices
 //-----------------------
 
 int new_up,new_right;
 int new_vertex;
 
 
 // Map that stores the set of rotations (as a pairs consisting of 
 // the up_equivalent and the right_equivalent) that move the vertex specified
 // by the first entry in key pair to the position of the second one:
 std::map<std::pair<int,int>, std::set<std::pair<int,int> > > 
  required_rotation;
 
 // Loop over all vertices
 for (int vertex=LDB;vertex<=RUF;vertex++)
  {
   
   new_up=U;
   new_right=R;
   new_vertex=OcTree::rotate(new_up,new_right,vertex);
   required_rotation[std::make_pair(vertex,new_vertex)].insert(
    std::make_pair(new_up,new_right));

    new_up=U;
    new_right=B;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=U;
    new_right=L;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=U;
    new_right=F;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    
    
   


    new_up=D;
    new_right=R;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=D;
    new_right=B;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=D;
    new_right=L;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=D;
    new_right=F;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));




    new_up=R;
    new_right=D;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=R;
    new_right=B;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=R;
    new_right=U;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=R;
    new_right=F;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));





    new_up=L;
    new_right=D;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=L;
    new_right=B;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=L;
    new_right=U;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=L;
    new_right=F;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));



    new_up=F;
    new_right=R;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=F;
    new_right=L;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=F;
    new_right=U;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=F;
    new_right=D;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));






    new_up=B;
    new_right=R;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=B;
    new_right=L;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=B;
    new_right=U;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));

    new_up=B;
    new_right=D;
    new_vertex=OcTree::rotate(new_up,new_right,vertex);
    required_rotation[std::make_pair(vertex,new_vertex)].insert(
     std::make_pair(new_up,new_right));


  }


 // Each vertex is part of three edges. This container stores the
 // vertices in each of the three edge neighbours that are 
 // adjacent to this node if there's no relative rotation between
 // the elements. 
 std::map<int,Vector<int> > vertex_in_edge_neighbour;
 
 
 // Each vertex is part of three edges. This container stores the
 // vertices at the other end of the edge
 std::map<int,Vector<int> > other_vertex_on_edge;
 
 // Each vertex is part of three edges. This container stores the
 // vertices in the adjacent element at the other end of the edge
 // assuming there are no rotations between the elements.
 std::map<int,Vector<int> > other_vertex_in_edge_neighbour;
 
 
 
         vertex_in_edge_neighbour[LDB].resize(3);
 vertex_in_edge_neighbour[LDB][0]=LUF; // Pattern: exactly one letter matches
 vertex_in_edge_neighbour[LDB][1]=RDF;
 vertex_in_edge_neighbour[LDB][2]=RUB;
 
         other_vertex_on_edge[LDB].resize(3);
 other_vertex_on_edge[LDB][0]=RDB; // Pattern: opposite of the matching letter
 other_vertex_on_edge[LDB][1]=LUB;
 other_vertex_on_edge[LDB][2]=LDF;
 
 other_vertex_in_edge_neighbour[LDB].resize(3);
 other_vertex_in_edge_neighbour[LDB][0]=RUF; // Pattern: full reflection
 other_vertex_in_edge_neighbour[LDB][1]=RUF;
 other_vertex_in_edge_neighbour[LDB][2]=RUF;
 
 
 
 
          vertex_in_edge_neighbour[RDB].resize(3);
  vertex_in_edge_neighbour[RDB][0]=RUF; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[RDB][1]=LDF;
  vertex_in_edge_neighbour[RDB][2]=LUB;

          other_vertex_on_edge[RDB].resize(3);
  other_vertex_on_edge[RDB][0]=LDB; // Pattern: opposite of the matching letter
  other_vertex_on_edge[RDB][1]=RUB;
  other_vertex_on_edge[RDB][2]=RDF;

          other_vertex_in_edge_neighbour[RDB].resize(3);
  other_vertex_in_edge_neighbour[RDB][0]=LUF; // Pattern: full reflection
  other_vertex_in_edge_neighbour[RDB][1]=LUF;
  other_vertex_in_edge_neighbour[RDB][2]=LUF;




          vertex_in_edge_neighbour[LUB].resize(3);
  vertex_in_edge_neighbour[LUB][0]=LDF; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[LUB][1]=RUF;
  vertex_in_edge_neighbour[LUB][2]=RDB;

          other_vertex_on_edge[LUB].resize(3);
  other_vertex_on_edge[LUB][0]=RUB; // Pattern: opposite of the matching letter
  other_vertex_on_edge[LUB][1]=LDB;
  other_vertex_on_edge[LUB][2]=LUF;

          other_vertex_in_edge_neighbour[LUB].resize(3);
  other_vertex_in_edge_neighbour[LUB][0]=RDF; // Pattern: full reflection
  other_vertex_in_edge_neighbour[LUB][1]=RDF;
  other_vertex_in_edge_neighbour[LUB][2]=RDF;




          vertex_in_edge_neighbour[RUB].resize(3);
  vertex_in_edge_neighbour[RUB][0]=RDF; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[RUB][1]=LUF;
  vertex_in_edge_neighbour[RUB][2]=LDB;

          other_vertex_on_edge[RUB].resize(3);
  other_vertex_on_edge[RUB][0]=LUB; // Pattern: opposite of the matching letter
  other_vertex_on_edge[RUB][1]=RDB;
  other_vertex_on_edge[RUB][2]=RUF;

          other_vertex_in_edge_neighbour[RUB].resize(3);
  other_vertex_in_edge_neighbour[RUB][0]=LDF; // Pattern: full reflection
  other_vertex_in_edge_neighbour[RUB][1]=LDF;
  other_vertex_in_edge_neighbour[RUB][2]=LDF;




          vertex_in_edge_neighbour[LDF].resize(3);
  vertex_in_edge_neighbour[LDF][0]=LUB; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[LDF][1]=RDB;
  vertex_in_edge_neighbour[LDF][2]=RUF;

          other_vertex_on_edge[LDF].resize(3);
  other_vertex_on_edge[LDF][0]=RDF; // Pattern: opposite of the matching letter
  other_vertex_on_edge[LDF][1]=LUF;
  other_vertex_on_edge[LDF][2]=LDB;

          other_vertex_in_edge_neighbour[LDF].resize(3);
  other_vertex_in_edge_neighbour[LDF][0]=RUB; // Pattern: full reflection
  other_vertex_in_edge_neighbour[LDF][1]=RUB;
  other_vertex_in_edge_neighbour[LDF][2]=RUB;




          vertex_in_edge_neighbour[RDF].resize(3);
  vertex_in_edge_neighbour[RDF][0]=RUB; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[RDF][1]=LDB;
  vertex_in_edge_neighbour[RDF][2]=LUF;

          other_vertex_on_edge[RDF].resize(3);
  other_vertex_on_edge[RDF][0]=LDF; // Pattern: opposite of the matching letter
  other_vertex_on_edge[RDF][1]=RUF;
  other_vertex_on_edge[RDF][2]=RDB;

          other_vertex_in_edge_neighbour[RDF].resize(3);
  other_vertex_in_edge_neighbour[RDF][0]=LUB; // Pattern: full reflection
  other_vertex_in_edge_neighbour[RDF][1]=LUB;
  other_vertex_in_edge_neighbour[RDF][2]=LUB;




          vertex_in_edge_neighbour[LUF].resize(3);
  vertex_in_edge_neighbour[LUF][0]=LDB; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[LUF][1]=RUB;
  vertex_in_edge_neighbour[LUF][2]=RDF;

          other_vertex_on_edge[LUF].resize(3);
  other_vertex_on_edge[LUF][0]=RUF; // Pattern: opposite of the matching letter
  other_vertex_on_edge[LUF][1]=LDF;
  other_vertex_on_edge[LUF][2]=LUB;

          other_vertex_in_edge_neighbour[LUF].resize(3);
  other_vertex_in_edge_neighbour[LUF][0]=RDB; // Pattern: full reflection
  other_vertex_in_edge_neighbour[LUF][1]=RDB;
  other_vertex_in_edge_neighbour[LUF][2]=RDB;




          vertex_in_edge_neighbour[RUF].resize(3);
  vertex_in_edge_neighbour[RUF][0]=RDB; // Pattern: exactly one letter matches
  vertex_in_edge_neighbour[RUF][1]=LUB;
  vertex_in_edge_neighbour[RUF][2]=LDF;

          other_vertex_on_edge[RUF].resize(3);
  other_vertex_on_edge[RUF][0]=LUF; // Pattern: opposite of the matching letter
  other_vertex_on_edge[RUF][1]=RDF;
  other_vertex_on_edge[RUF][2]=RUB;

          other_vertex_in_edge_neighbour[RUF].resize(3);
  other_vertex_in_edge_neighbour[RUF][0]=LDB; // Pattern: full reflection
  other_vertex_in_edge_neighbour[RUF][1]=LDB;
  other_vertex_in_edge_neighbour[RUF][2]=LDB;




  // Loop over all vertices in the reference element
  for (int vertex=LDB;vertex<=RUF;vertex++)
   {

    // Loop over the three edges that are connected to this vertex
    for (unsigned i=0;i<3;i++)
     {
      
      // This is the other vertex along this edge
      int other_vertex=other_vertex_on_edge[vertex][i];

      // This is the vertex in the edge neighbour that corresponds
      // to the present vertex in the reference element (in
      // the absence of rotations)
      int unrotated_neigh_vertex=vertex_in_edge_neighbour[vertex][i];

      // This is the vertex in the edge neighbour that corresponds
      // to the other vertex in the reference element (in
      // the absence of rotations)
      int unrotated_neigh_other_vertex=
       other_vertex_in_edge_neighbour[vertex][i];

      // Loop over all vertices in the neighbour element
      for (int neigh_vertex=LDB;neigh_vertex<=RUF;neigh_vertex++)
       {
        
        // What rotations would turn the neigh_vertex
        // into the unrotated_neigh_vertex?
        std::set<std::pair<int,int> > vertex_rot=
         required_rotation[std::make_pair(neigh_vertex,
                                          unrotated_neigh_vertex)];
                   


        // Loop over all "other" vertices in the neighbour element
        for (int neigh_other_vertex=LDB;neigh_other_vertex<=RUF;
             neigh_other_vertex++)
         {
                    
          // What rotations would turn the other_neigh_vertex
          // into the unrotated_other_neigh_vertex?
          std::set<std::pair<int,int> > other_vertex_rot=
           required_rotation[std::make_pair(neigh_other_vertex,
                                            unrotated_neigh_other_vertex)];
          
          // What are the common rotations?
          std::set<std::pair<int,int> > common_rotations;
          
          // Get the intersection of the two sets
          std::set_intersection(vertex_rot.begin(),
                                vertex_rot.end(),
                                other_vertex_rot.begin(),
                                other_vertex_rot.end(),
                                inserter(common_rotations,
                                         common_rotations.begin()));
          

          if (common_rotations.size()>0)
           {
            for (std::set<std::pair<int,int> >::iterator it=
                  common_rotations.begin();
                 it!=common_rotations.end();it++)
             {
              // Copy into container

              // First: up equivalent:
              Up_and_right_equivalent_for_pairs_of_vertices[
               std::make_pair(std::make_pair(vertex,neigh_vertex),
                              std::make_pair(other_vertex,
                                             neigh_other_vertex))].
               first=it->first;
              
              // Second: Right equivalent
              Up_and_right_equivalent_for_pairs_of_vertices[
               std::make_pair(std::make_pair(vertex,neigh_vertex),
                              std::make_pair(other_vertex,
                                             neigh_other_vertex))].
               second=it->second;

             }
           }          
          
         }
       }
     }
   }
  


}



//================================================================
/// \short Is the edge neighbour (for edge "edge")  specified via the pointer 
/// also a face neighbour for one of the two adjacent faces?
//================================================================
bool OcTree::edge_neighbour_is_face_neighbour(const int& edge,
                                              OcTree* edge_neigh_pt) const
{

#ifdef PARANOID
 // No paranoid check needed -- the default for the switch statement 
 // catches illegal values for edge
#endif


 // Catch stupid case: Null doesn't have a face neighbour...
 if (edge_neigh_pt==0) 
  {return false;}

 using namespace OcTreeNames;

 // Auxiliary variables
 int face;
 Vector<unsigned> translate_s(3);
 Vector<double> s_sw(3);
 Vector<double> s_ne(3);
 int reflected_face;
 int diff_level;

 OcTree* face_neigh_pt=0;

 switch (edge)
  {
   case LB:

    // Get first face neighbour
    face=L;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if(face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=B;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
    
  case RB:


    // Get first face neighbour
    face=R;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=B;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
    
  case DB:

    // Get first face neighbour
    face=D;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=B;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    

  case UB:

    // Get first face neighbour
    face=U;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=B;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
  case LD:
  

    // Get first face neighbour
    face=L;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=D;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
  case RD: 
   

    // Get first face neighbour
    face=R;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=D;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
  case LU: 

    // Get first face neighbour
    face=L;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=U;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    

  case RU:

    // Get first face neighbour
    face=R;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
      
     }

    // Get second face neighbour
    face=U;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
    
  case LF:
  

    // Get first face neighbour
    face=L;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     { 
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }
      
    // Get second face neighbour
    face=F;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
  case RF:
  
    // Get first face neighbour
    face=R;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    // Get second face neighbour
    face=F;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    

  case DF:
  
    // Get first face neighbour
    face=D;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }


    // Get second face neighbour
    face=F;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
    
    
  case UF:

    // Get first face neighbour
    face=U;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }


    // Get second face neighbour
    face=F;
    face_neigh_pt=gteq_face_neighbour(face,
                                      translate_s,
                                      s_sw,
                                      s_ne, 
                                      reflected_face,
                                      diff_level);
    // Check if they agree...
    if (face_neigh_pt!=0)
     {
      if (face_neigh_pt==edge_neigh_pt)
       {
        return true;
       }
     }

    break;
   
  default:

   // There is no face neighbour in this direction so they can't 
   // agree:
   std::ostringstream error_stream;
   error_stream << "Never get here! Edge:" << Direct_string[edge]
                << " " << edge << std::endl;
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);


  }

 // If we've made it to here then we've located the requested edge
 // but found that none of its two face neighbours match the specified
 // edge neighbour:
 return false;
 

}


//================================================================
/// Find (pointer to) `greater-or-equal-sized face neighbour' in
/// given direction (L/R/U/D/F/B).
/// Another way of interpreting this is that we're looking for
/// the neighbour across the present element's face 'direction'.
/// The various arguments return additional information about the
/// size and relative orientation of the neighbouring octree.
/// To interpret these we use the following
/// <B>General convention:</B>
/// - Each face of the element that is represented by the octree
///   is parametrised by two (of the three)
///   local coordinates that parametrise the entire 3D element. E.g.
///   the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
///   is parametrised by (s[0],s[2]); etc. We always identify the
///   in-face coordinate with the lower (3D) index with the subscript
///   _lo and the one with the larger (3D) index with the subscript _hi.
/// .
/// With this convention, the interpretation of the arguments is
/// as follows:
/// - The vector \c translate_s turns the index of the local coordinate
///   in the present octree into that of the neighbour. If there are no
///   rotations then \c translate_s[i] = i.
/// - In the present octree, the "south west" vertex of the face
///   between the present octree and its neighbour is located at
///   S_lo=-1, S_hi=-1. This point is located at the (3D) local
///   coordinates (\c s_sw[0], \c s_sw[1], \c s_sw[2])
///   in the neighbouring octree.
/// - ditto with s_ne: In the present octree, the "north east" vertex
///   of the face between the present octree and its neighbour is located at
///   S_lo=+1, S_hi=+1. This point is located
///   at the (3D) local coordinates (\c s_ne[0], \c s_ne[1], \c s_ne[2])
///   in the neighbouring octree.
/// - We're looking for a neighbour in the specified \c direction. When
///   viewed from the neighbouring octree, the face that separates
///   the present octree from its neighbour is the neighbour's face
///   \c face. If there's no rotation between the two octrees, this is a
///   simple reflection: For instance, if we're looking
///   for a neighhbour in the \c R [ight] \c direction, \c face will
///   be \c L [eft]
/// - \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 octree.
//=================================================================
OcTree* OcTree::gteq_face_neighbour(const int& direction,
                                    Vector<unsigned> &translate_s,
                                    Vector<double>& s_sw,
                                    Vector<double>& s_ne,
                                    int& face,
                                    int& diff_level) const
{

 using namespace OcTreeNames;

#ifdef PARANOID
    if ((direction!=L)&&(direction!=R)&&
        (direction!=F)&&(direction!=B)&&
        (direction!=U)&&(direction!=D))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong direction: " << Direct_string[direction] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif


 // 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:
 OcTreeRoot* orig_root_pt = dynamic_cast<OcTreeRoot*>(Root_pt);

 // Initialise offset in local coordinate
 double s_difflo=0;
 double s_diffhi=0;
 
 // Initialise difference in level
 diff_level=0;

 // Find neighbour
 OcTree* return_pt=gteq_face_neighbour(direction,s_difflo,s_diffhi,
                                       diff_level,max_level,
                                       orig_root_pt);

 OcTree* neighb_pt=return_pt;

 //Initialise the translation scheme
 for(unsigned i=0;i<3;i++) {translate_s[i] = i;}

 // If neighbour exists: What's the direction of the interfacial
 // face when viewed from within the neighbour element?
 if (neighb_pt!=0)
  {
   //Find the reflection of the original direction, which will be the 
   //direction to the face in the neighbour, if there are no rotations.
   int reflected_dir = Reflect_face[direction];

   // These coordinates are the coordinates of the ne and sw points
   // in the neighbour (provided there are no rotations -- we'll correct
   // for these below)
   s_sw[0]=S_base(0,reflected_dir)+
    S_steplo(0,reflected_dir)*s_difflo+S_stephi(0,reflected_dir)*s_diffhi;
   
   s_sw[1]=S_base(1,reflected_dir)+
    S_steplo(1,reflected_dir)*s_difflo+S_stephi(1,reflected_dir)*s_diffhi;
   
   s_sw[2]=S_base(2,reflected_dir)+
    S_steplo(2,reflected_dir)*s_difflo+S_stephi(2,reflected_dir)*s_diffhi;
   
   s_ne[0]=S_base(0,reflected_dir)+
    S_steplo(0,reflected_dir)*pow(2.0,diff_level)+
    S_steplo(0,reflected_dir)*s_difflo+
    S_stephi(0,reflected_dir)*pow(2.0,diff_level)+
    S_stephi(0,reflected_dir)*s_diffhi;
   
   s_ne[1]=S_base(1,reflected_dir)+
    S_steplo(1,reflected_dir)*pow(2.0,diff_level)
    +S_steplo(1,reflected_dir)*s_difflo+
    S_stephi(1,reflected_dir)*pow(2.0,diff_level)
    +S_stephi(1,reflected_dir)*s_diffhi;
   
   s_ne[2]=S_base(2,reflected_dir)+
    S_steplo(2,reflected_dir)*pow(2.0,diff_level)
    +S_steplo(2,reflected_dir)*s_difflo+
    S_stephi(2,reflected_dir)*pow(2.0,diff_level)
    +S_stephi(2,reflected_dir)*s_diffhi;
   
   //If there is no rotation then the new direction is the same as the
   //old direction
   int new_dir = direction;
   
   // If necessary, rotate things around (the orientation of Up in the 
   // neighbour might be be different from that in the present element)
   // If the root of the neighbour is the not same as the one of the present
   // element then their orientations may not be the same and the new direction
   // is given by :
   if (neighb_pt->Root_pt!=Root_pt)
    {
     new_dir=rotate(orig_root_pt->up_equivalent(
                     neighb_pt->Root_pt),
                    orig_root_pt->right_equivalent(
                     neighb_pt->Root_pt),
                    direction);
    }

   // What's the direction of the interfacial face when viewed from within
   // the neighbour element?
   face=Reflect_face[new_dir];
   
   Vector<double> s_sw_new(3), s_ne_new(3);
   
   // If the root of the present element is different from the root 
   // of his neighbour, we have to rotate the ne and sw coordinates
   // to have their equivalents in the neighbour's point of view.
   
    if (neighb_pt->Root_pt!=Root_pt)
     {
      int tmp_dir;     
      Vector<int> vect1(3);
      Vector<int> vect2(3);
      Vector<int> vect3(3);
      DenseMatrix<int> Mat_rot(3);
      
      
      // All this is just to compute the rotation matrix 
      tmp_dir=rotate(orig_root_pt->
                      up_equivalent(neighb_pt->Root_pt),
                      orig_root_pt->
                      right_equivalent(neighb_pt->Root_pt),
                      R);
      vect1=Direction_to_vector[tmp_dir];
      
      
      tmp_dir=rotate(orig_root_pt->
                      up_equivalent(neighb_pt->Root_pt),
                      orig_root_pt->
                      right_equivalent(neighb_pt->Root_pt),
                      U);
      vect2=Direction_to_vector[tmp_dir];
      
      
      tmp_dir=rotate(orig_root_pt->
                      up_equivalent(neighb_pt->Root_pt),
                      orig_root_pt->
                      right_equivalent(neighb_pt->Root_pt),
                      F);
      vect3=Direction_to_vector[tmp_dir];
      
      
      // Setup the inverse rotation matrix
      for(int i=0; i<3; i++)
       {
        Mat_rot(i,0)=vect1[i];
        Mat_rot(i,1)=vect2[i];
        Mat_rot(i,2)=vect3[i];
       }
      
      
      //Initialise the translation scheme
      Vector<int> translate_s_new(3);
      
      // Then the rotation of the coordinates
      for(unsigned i=0;i<3;i++)
       {
        s_ne_new[i] = 0.0;
        s_sw_new[i] = 0.0;
        translate_s_new[i] = 0;
        for(unsigned k=0;k<3;k++) 
         {
          s_ne_new[i] += Mat_rot(i,k)*s_ne[k];
          s_sw_new[i] += Mat_rot(i,k)*s_sw[k]; 
          translate_s_new[i] += Mat_rot(i,k)*translate_s[k];
         }
       }
      
      s_ne=s_ne_new;
      s_sw=s_sw_new;
      
      //Set the translation scheme
      for(unsigned i=0;i<3;i++)
       {
        // abs is ok here; not fabs!
        translate_s[i] = std::abs(translate_s_new[i]);
       }
     }
    
  } // end of if(neighb_pt!=0)
 
 return return_pt;
 
}



//================================================================
///  Find (pointer to) `greater-or-equal-sized true edge neighbour' in 
/// the given direction (LB,RB,DB,UB [the back edges],
/// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]). 
/// Another way of interpreting this is that we're looking for 
/// the neighbour across the present element's edge 'direction'.
/// The various arguments return additional information about the
/// size and relative orientation of the neighbouring octree.
/// Each edge of the element that is represented by the octree 
/// is parametrised by one (of the three) local coordinates that 
/// parametrise the entire 3D element. E.g. the L[eft]B[ack] edge 
/// is parametrised by s[1]; the "low" vertex of this edge 
/// (located at the low value of this coordinate, i.e. at s[1]=-1)
/// is L[eft]D[own]B[ack]. The "high" vertex of this edge (located
/// at the high value of this coordinate, i.e. at s[1]=1) is
/// L[eft]U[p]B[ack]; etc
/// The interpretation of the arguments is as follows:
/// - In a forest, an OcTree can have multiple edge neighbours
///   (across an edge where multiple trees meet). \c i_root_edge_neighbour
///   specifies which of these is used. Use this as "reverse communication":
///   First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
///   initialised to anything you want (zero, ideally). On return from
///   the fct, \c n_root_edge_neighour contains the total number of true
///   edge neighbours, so additional calls to the fct with 
///   \c i_root_edge_neighbour>0 can be made until they've all been visited.
/// - The vector \c translate_s turns the index of the local coordinate
///   in the present octree into that of the neighbour. If there are no
///   rotations then \c translate_s[i] = i.
/// - The "low" vertex of the edge in the present octree
///   coincides with a certain vertex in the edge neighbour.
///   In terms of the neighbour's local coordinates, this point is 
///   located at the (3D) local coordinates (\c s_lo[0], \c s_lo[1], 
///   \c s_lo[2])
/// - ditto with s_hi: The "high" vertex of the edge in the present octree
///   coincides with a certain vertex in the edge neighbour.
///   In terms of the neighbour's local coordinates, this point is 
///   located at the (3D) local coordinates (\c s_hi[0], \c s_hi[1], 
///   \c s_hi[2])
/// - We're looking for a neighbour in the specified \c direction. When
///   viewed from the neighbouring octree, the edge that separates 
///   the present octree from its neighbour is the neighbour's edge 
///   \c edge. If there's no rotation between the two octrees, this is a 
///   simple reflection: For instance, if we're looking
///   for a neighhbour in the \c DB \c direction, \c edge will 
///   be \c UF.
/// - \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 octree.
/// .
/// \b Important: We're only looking for \b true edge neighbours
/// i.e. edge neigbours that are not also face neighbours. This is an
/// important difference to Samet's terminology. If the neighbour
/// in a certain direction is not a true edge neighbour, or if there
/// is no neighbour, then this function returns NULL.
//=================================================================
OcTree* OcTree::gteq_true_edge_neighbour(const int& direction,
                                         const unsigned& i_root_edge_neighbour,
                                         unsigned& nroot_edge_neighbour,
                                         Vector<unsigned> &translate_s,
                                         Vector<double>& s_lo,
                                         Vector<double>& s_hi,
                                         int& edge,
                                         int& diff_level) const
{
 
 using namespace OcTreeNames;


#ifdef PARANOID
    if ((direction!=LB)&&(direction!=RB)&&(direction!=DB)&&(direction!=UB)&&
        (direction!=LD)&&(direction!=RD)&&(direction!=LU)&&(direction!=RU)&&
        (direction!=LF)&&(direction!=RF)&&(direction!=DF)&&(direction!=UF))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong direction: " << Direct_string[direction] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

 // 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:
 OcTreeRoot* orig_root_pt = dynamic_cast<OcTreeRoot*>(Root_pt);

 // Initialise offset in local coordinate along edge
 double s_diff=0;
 
 // Initialise difference in level
 diff_level=0;

 // Find edge neighbour
 OcTree* return_pt=gteq_edge_neighbour(direction,
                                       i_root_edge_neighbour,
                                       nroot_edge_neighbour,
                                       s_diff,
                                       diff_level,max_level,
                                       orig_root_pt);

 // Only use "true" edge neighbours
 if (edge_neighbour_is_face_neighbour(direction,return_pt))
  {
   return_pt=0;
  }

 // By default, we return what was returned as the true edge neighbour.
 OcTree* neighb_pt=return_pt;
 
 //Initialise the translation scheme
 for(unsigned i=0;i<3;i++) {translate_s[i] = i;}

 // If neighbour exists: What's the direction of the interfacial
 // edge when viewed from within the neighbour element?
 if (neighb_pt!=0)
  {
   //Find the reflection of the original direction, which will be the 
   //direction to the edge in the neighbour, if there are no rotations.
   int reflected_dir = Reflect_edge[direction];

   // These coordinates are the coordinates of the "low" and "high" points 
   // in the neighbour (provided there are no rotations -- we'll correct
   // for these below)
   s_lo[0]=S_base_edge(0,reflected_dir)+
    S_step_edge(0,reflected_dir)*s_diff;
   
   s_lo[1]=S_base_edge(1,reflected_dir)+
    S_step_edge(1,reflected_dir)*s_diff;
   
   s_lo[2]=S_base_edge(2,reflected_dir)+
    S_step_edge(2,reflected_dir)*s_diff;
   
   s_hi[0]=S_base_edge(0,reflected_dir)+
    S_step_edge(0,reflected_dir)*pow(2.0,diff_level)+
    S_step_edge(0,reflected_dir)*s_diff;
   
   s_hi[1]=S_base_edge(1,reflected_dir)+
    S_step_edge(1,reflected_dir)*pow(2.0,diff_level)
    +S_step_edge(1,reflected_dir)*s_diff;
   
   s_hi[2]=S_base_edge(2,reflected_dir)+
    S_step_edge(2,reflected_dir)*pow(2.0,diff_level)
    +S_step_edge(2,reflected_dir)*s_diff;
   
   //If there is no rotation then the new direction is the same as the
   //old direction
   int new_dir = direction;
   

   // If necessary, rotate things around (the orientation of Up in the 
   // neighbour might be be different from that in the present element)
   // If the root of the neighbour is the not same as the one of the present
   // element then their orientations may not be the same and the new 
   // direction is given by :
   if (neighb_pt->Root_pt!=Root_pt)
    {     
     int new_up=orig_root_pt->
      up_equivalent(neighb_pt->Root_pt);

     int new_right=orig_root_pt->
      right_equivalent(neighb_pt->Root_pt);
     
     new_dir=rotate(new_up,new_right,direction);
    }
      
   // What's the direction of the interfacial edge when viewed from within
   // the neighbour element (including rotations!)
   edge=Reflect_edge[new_dir];
   
   // Get ready to rotate the local coordinates
   Vector<double> s_lo_new(3), s_hi_new(3);
   
   // If the root of the present element is different from the root 
   // of his neighbour, we have to rotate the lo and hi coordinates
   // to have their equivalents from the neighbour's point of view.
   if ((neighb_pt->Root_pt!=Root_pt)) 
    {

     int tmp_dir;     
     Vector<int> vect1(3);
     Vector<int> vect2(3);
     Vector<int> vect3(3);
     DenseMatrix<int> Mat_rot(3);
          
     // All this is just to compute the rotation matrix 
     tmp_dir = rotate(orig_root_pt->up_equivalent(
                               neighb_pt->Root_pt),
                              orig_root_pt->right_equivalent(
                               neighb_pt->Root_pt), 
                              R);
     vect1=Direction_to_vector[tmp_dir];
     
     
     tmp_dir = rotate(orig_root_pt->up_equivalent(
                               neighb_pt->Root_pt),
                              orig_root_pt->right_equivalent(
                               neighb_pt->Root_pt),
                              U);
     vect2=Direction_to_vector[tmp_dir];
     
     
     tmp_dir = rotate(orig_root_pt->up_equivalent(
                               neighb_pt->Root_pt),
                              orig_root_pt->right_equivalent(
                               neighb_pt->Root_pt),
                              F);
     vect3=Direction_to_vector[tmp_dir];
     

     // Setup the inverse rotation matrix
     for(int i=0; i<3; i++)
      {
       Mat_rot(i,0)=vect1[i];
       Mat_rot(i,1)=vect2[i];
       Mat_rot(i,2)=vect3[i];
      }
     
     //Initialise the translation scheme
     Vector<int> translate_s_new(3);
     
     // Then the rotation of the coordinates
     for(unsigned i=0;i<3;i++)
      {
       s_hi_new[i] = 0.0;
       s_lo_new[i] = 0.0;
       translate_s_new[i] = 0;
       for(unsigned k=0;k<3;k++) 
        {
         s_hi_new[i] += Mat_rot(i,k)*s_hi[k];
         s_lo_new[i] += Mat_rot(i,k)*s_lo[k]; 
         translate_s_new[i] += Mat_rot(i,k)*translate_s[k];
        }
      }
     
     s_lo=s_lo_new;
     s_hi=s_hi_new;
     
     //Set the translation scheme
     for(unsigned i=0;i<3;i++)
      {
       // abs is ok here; not fabs!
       translate_s[i] = std::abs(translate_s_new[i]);
      }
    }
   
  } // end if for (neighb_pt!=0)

 return return_pt;
 
}




//==========================================================================
/// Find `greater-or-equal-sized face neighbour' in given direction
/// (L/R/U/D/B/F).
///
/// This is an auxiliary routine which allows neighbour finding in adjacent
/// octrees. Needs to keep track of previous son types and
/// the maximum level to which search is performed.
///
/// Parameters:
///
/// - direction: (L/R/U/D/B/F) Direction in which neighbour has to be found.
/// - s_difflo/s_diffhi: Offset of left/down/back 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.
/// - face: We're looking for the neighbour across our face 'direction'
///   (L/R/U/D/B/F). When viewed from the neighbour, this face is
///   `face' (L/R/U/D/B/F). [If there's no relative rotation between neighbours
///    then this is a mere reflection, e.g. direction=F  --> face=B etc.]
/// - diff_level <= 0 indicates the difference in octree 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 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.
//===========================================================================
OcTree* OcTree::gteq_face_neighbour(const int& direction,
                                    double& s_difflo,
                                    double& s_diffhi,
                                    int& diff_level,
                                    int max_level,
                                    OcTreeRoot* orig_root_pt) const
{

using namespace OcTreeNames;


#ifdef PARANOID
    if ((direction!=L)&&(direction!=R)&&
        (direction!=F)&&(direction!=B)&&
        (direction!=U)&&(direction!=D))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong direction: " << Direct_string[direction] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

 OcTree* next_el_pt;
 OcTree* return_el_pt;


 // STEP 1: Find the neighbour's father
 //--------

 // Does the element have a father?
 if (Father_pt!=0)
  {

   // If the present octant (whose location inside its
   // father element is specified by Son_type) is adjacent to the
   // father's face 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
   if (Is_adjacent(direction,Son_type))
    {

     next_el_pt=dynamic_cast<OcTree*>(Father_pt)->
      gteq_face_neighbour(direction,
                          s_difflo, s_diffhi,diff_level,max_level,
                          orig_root_pt);
    }
   
   
   // If the present octant is not adjacent to the
   // father's face in the required direction, then the
   // neighbour has the same father and will be obtained
   // by the appropriate reflection inside the father element
   else
    {
     next_el_pt = dynamic_cast<OcTree*>(Father_pt);
    }


   // We're about to ascend one level:
   diff_level-=1;

   // Work out position of lower-left corner of present face
   // in its father element

   s_difflo+=pow(0.5,-diff_level)*S_directlo(direction,Son_type);
   s_diffhi+=pow(0.5,-diff_level)*S_directhi(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.
     else
       {
         
        // By default (in the absence of rotations) we obtain the
        // son octant by reflecting
        int son_octant=Reflect(direction,Son_type);
         

        // If there is a rotation, we rotate the son octant
        if (orig_root_pt!=next_el_pt->Root_pt)
           {
            int my_up=dynamic_cast<OcTreeRoot*>(Root_pt)
             ->up_equivalent(next_el_pt->Root_pt);

            int my_right=dynamic_cast<OcTreeRoot*>(Root_pt)
             ->right_equivalent(next_el_pt->Root_pt);

            son_octant=rotate(my_up,my_right,son_octant);
           }
        
        return_el_pt=dynamic_cast<OcTree*>(next_el_pt->Son_pt[son_octant]);
        
        
        // Work out position of lower-left corner of present face
        // in next higher element
        s_difflo-=pow(0.5,-diff_level)*S_directlo(direction,Son_type);
        s_diffhi-=pow(0.5,-diff_level)*S_directhi(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)
    {
     return_el_pt=dynamic_cast<OcTree*>(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;
 
}




//==========================================================================
/// Find `greater-or-equal-sized edge neighbour' in given direction
///  (LB,RB,DB,UB [the back edges],
/// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]). 
///
/// This is an auxiliary routine which allows neighbour finding in adjacent
/// octrees. Needs to keep track of previous son types and
/// the maximum level to which search is performed.
///
/// Parameters:
///
/// - direction: (LB,RB/...) Direction in which neighbour has to be found.
/// - In a forest, an OcTree can have multiple edge neighbours
///   (across an edge where multiple trees meet). \c i_root_edge_neighbour
///   specifies which of these is used. Use this as "reverse communication":
///   First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
///   initialised to anything you want (zero, ideally). On return from
///   the fct, \c n_root_edge_neighour contains the total number of true
///   edge neighbours, so additional calls to the fct with 
///   \c i_root_edge_neighbour>0 can be made until they've all been visited.
/// - s_diff: Offset of left/down/back 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.
/// - diff_level <= 0 indicates the difference in octree 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 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.
//===========================================================================
OcTree* OcTree::gteq_edge_neighbour(const int& direction,
                                    const unsigned& i_root_edge_neighbour,
                                    unsigned& nroot_edge_neighbour,
                                    double& s_diff,
                                    int& diff_level,
                                    int max_level,
                                    OcTreeRoot* orig_root_pt) const
{

using namespace OcTreeNames;


#ifdef PARANOID
 if ((direction!=LB)&&(direction!=RB)&&(direction!=DB)&&(direction!=UB)&&
     (direction!=LD)&&(direction!=RD)&&(direction!=LU)&&(direction!=RU)&&
     (direction!=LF)&&(direction!=RF)&&(direction!=DF)&&(direction!=UF))
  {
   std::ostringstream error_stream;
   error_stream << "Wrong direction: " << Direct_string[direction] 
                << std::endl;
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Initialise total number of edge neighbours available across
 // edges of octree roots
 nroot_edge_neighbour=0;


 OcTree* next_el_pt=0;
 OcTree* return_el_pt=0;


 // STEP 1: Find the common ancestor
 //--------

 // Does the element have a father?
 if (Father_pt!=0)
  {


   // If the present octant (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
   if (Is_adjacent(direction,Son_type))
    {
     next_el_pt=dynamic_cast<OcTree*>(Father_pt)->
      gteq_edge_neighbour(direction,
                          i_root_edge_neighbour,
                          nroot_edge_neighbour,
                          s_diff,diff_level,max_level,
                          orig_root_pt);
    }
   // If the present octant (whose location inside its
   // father element is specified by Son_type) is adjacent to the
   // father's face 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,
   // crossing the face as we do so.
   else if (Common_face(direction,Son_type)!=OMEGA)
    {
     // We're going across a face:
     double s_difflo=0.0;
     double s_diffhi=0.0;
     int diff_level_edge=0;

     next_el_pt=dynamic_cast<OcTree*>(Father_pt)->
      gteq_face_neighbour(Common_face(direction,Son_type),
                          s_difflo, s_diffhi,
                          diff_level_edge,max_level,
                          orig_root_pt);
    }

   // If the present octant is not adjacent to the
   // father's face/edge in the required direction, then the
   // neighbour has the same father and will be obtained
   // by the appropriate reflection inside the father element
   else
    {
     next_el_pt = dynamic_cast<OcTree*>(Father_pt);     
    }
   
   
   // We're about to ascend one level:
   diff_level-=1;
   
   // Work out position of "low" vertex of present edge
   // in its father element
   s_diff+=pow(0.5,-diff_level)*S_direct_edge(direction,Son_type);
   
   // STEP 2:  We have now located the neighbour's father and need to
   // -------  find the appropriate son.
   
   // Buffer cases where ...
   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.
     else
      {
       
       // By default (in the absence of rotations) we obtain the
       // son octant by reflecting
       int son_octant=Reflect(direction,Son_type);
       
       // If there is a rotation, we rotate the son octant
       if (orig_root_pt!=next_el_pt->Root_pt)
        {
         int my_up=dynamic_cast<OcTreeRoot*>(Root_pt)
          ->up_equivalent(next_el_pt->Root_pt);
         
         int my_right=dynamic_cast<OcTreeRoot*>(Root_pt)
          ->right_equivalent(next_el_pt->Root_pt);
         
         son_octant=rotate(my_up,my_right,son_octant);
        }
              
       return_el_pt=dynamic_cast<OcTree*>(next_el_pt->Son_pt[son_octant]);
              
       // Work out position of "low" vertex of present edge
       // in next higher element
       s_diff-=pow(0.5,-diff_level)*S_direct_edge(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
   {
    // Get total number of edge neighbours available across
    // edges of octree roots for return
    nroot_edge_neighbour=
     dynamic_cast<OcTreeRoot*>(Root_pt)->nedge_neighbour(direction);
    
    // Get vector of edge neighbours (if any) in appropriate direction 
    Vector<TreeRoot*> edge_neighbour_pt=
     dynamic_cast<OcTreeRoot*>(Root_pt)->edge_neighbour_pt(direction);
    unsigned n_neigh=edge_neighbour_pt.size();
    if (n_neigh>0)
     {
      return_el_pt=dynamic_cast<OcTree*>(
       edge_neighbour_pt[i_root_edge_neighbour]);
     }
    else
     {
      return_el_pt=0;
     }      
   }
 
 return return_el_pt;
 
}




//================================================================
/// 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 OcTree::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 vertices and their opposite points
 std::ofstream neighbours_file;
 std::ofstream no_true_edge_file;
 std::ofstream neighbours_txt_file;
 
 double max_error_face=0.0;
 OcTree::doc_face_neighbours(all_nodes_pt,neighbours_file,
                             neighbours_txt_file, max_error_face);
 
 double max_error_edge=0.0;
 OcTree::doc_true_edge_neighbours(all_nodes_pt,neighbours_file,
                                  no_true_edge_file,
                                  neighbours_txt_file, max_error_edge);
 

 bool failed=false;
 if (max_error_face>Max_neighbour_finding_tolerance)
  {
   oomph_info 
    << "\n \n Failed self_test() for OcTree because of faces: Max. error "
    << max_error_face << std::endl<< std::endl;
   failed=true;
  }

 if (max_error_edge>Max_neighbour_finding_tolerance)
  {
   oomph_info 
    << "\n \n Failed self_test() for OcTree because of edges: Max. error "
    << max_error_edge << std::endl<< std::endl;
   failed=true;
  }

 double max_error=max_error_face;
 if (max_error_edge>max_error) max_error=max_error_edge;

 if (failed)
  {
   return 1;
  }
 else
  {
   oomph_info << "Passed self_test() for OcTree: Max. error "
              << max_error << std::endl;
   return 0;
  }
 
}




//=================================================================
/// Doc/check all face neighbours of octree (nodes) contained in the
/// Vector forest_node_pt. Output into neighbours_file which can
/// be viewed from tecplot with OcTreeNeighbours.mcr
/// Neighbour info and errors are displayed on 
/// neighbours_txt_file.  Finally, compute the max. error between 
/// vertices when viewed from neighhbouring element. 
/// If the two filestreams are closed, output is suppressed. 
/// (Static function.)
//=================================================================
void OcTree::doc_face_neighbours(Vector<Tree*> forest_nodes_pt,
                                 std::ofstream& neighbours_file,
                                 std::ofstream& neighbours_txt_file,
                                 double& max_error)
{
 using namespace OcTreeNames;

 int diff_level;
 int face=OMEGA;
 
 Vector<double> s(3);
 Vector<double> x(3);

 Vector<double> s_sw(3);
 Vector<double> s_ne(3);
 Vector<unsigned> translate_s(3);

 Vector<double> x_small(3);
 Vector<double> x_large(3);

 
 // Initialise error in vertex positions
 max_error=0.0;

 // Loop over all elements
 // ----------------------
 unsigned long num_nodes=forest_nodes_pt.size();

 for (unsigned long i=0;i<num_nodes;i++)
  {
    

   // Doc the element itself
   OcTree* el_pt = dynamic_cast<OcTree*>(forest_nodes_pt[i]);

   //If the object is incomplete omit
   if(el_pt->object_pt()->nodes_built())
    {
     // Print it
     if (neighbours_file.is_open())
      {
       neighbours_file << "#---------------------------------" << std::endl;
       neighbours_file << "#The element itself: " << i << std::endl;
       neighbours_file << "#---------------------------------" << std::endl;
       dynamic_cast<RefineableQElement<3>*>(
        el_pt->object_pt())->output_corners(neighbours_file,"BLACK");
      }

   // Loop over directions to find neighbours
   // ----------------------------------------
   for (int direction=L;direction<=F;direction++)
    {
      
      // Initialise difference in levels and coordinate offset
     diff_level=0;

     // Find greater-or-equal-sized neighbour...
      OcTree* neighb_pt=el_pt->gteq_face_neighbour(direction, translate_s,
                                                   s_sw, s_ne,
                                                   face,diff_level);

      // If neighbour exists and nodes are created
      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
                            << ". Inside the neighbour the face is "
                            << Direct_string[face] <<  std::endl;
        }

       // Plot neighbour in the appropriate colour
       if (neighbours_file.is_open())
        {
         neighbours_file << "#---------------------------------" << std::endl;
         neighbours_file << "#Neighbour element: " << Direct_string[direction]
                         << std::endl;
         neighbours_file << "#---------------------------------" << std::endl;
         dynamic_cast<RefineableQElement<3>*>(
          neighb_pt->object_pt())->
          output_corners(neighbours_file,Colour[direction]);
        }

       // Check that local coordinates in the larger element
       // lead to the same spatial point as the node vertices
       // in the current element
       if (neighbours_file.is_open())
        {
         neighbours_file << "ZONE I=2  C=" << Colour[direction] << std::endl;
        }

       // "South west" vertex in the interfacial face
       //--------------------------------------------

       // Get coordinates in large (neighbour) element
       s[0]=s_sw[0];
       s[1]=s_sw[1];
       s[2]=s_sw[2];
       neighb_pt->object_pt()->get_x(s,x_large);

       // Get coordinates in small element
       Vector<double> s(3);
       s[0]=S_base(0,direction);
       s[1]=S_base(1,direction);
       s[2]=S_base(2,direction);
       el_pt->object_pt()->get_x(s,x_small);

       double error=
        sqrt(pow(x_small[0]-x_large[0],2)+
             pow(x_small[1]-x_large[1],2)+
             pow(x_small[2]-x_large[2],2));

       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);
        }


       // Check error and doc mismatch if required
       bool stop=false;
       std::ofstream mismatch_file;
       if (std::fabs(error)>Max_neighbour_finding_tolerance)
        {
         stop=true;
         mismatch_file.open("mismatch.dat");
         mismatch_file << "ZONE" << std::endl;
         mismatch_file<< x_large[0]<<" "<<x_large[1]<<" "<<x_large[2]<<" 2 \n";
         mismatch_file<< x_small[0]<<" "<<x_small[1]<<" "<<x_small[2]<<" 3 \n";
        }

       if (neighbours_file.is_open())
        {
         neighbours_file << "#SOUTH WEST: " << std::endl;
         neighbours_file << x_large[0] << " "
                         << x_large[1] << " "
                         << x_large[2] <<"  40 \n";
        }


       // "North east" vertex in the interfacial face
       //--------------------------------------------

       // Get coordinates in large (neighbour) element
       s[0]=s_ne[0];
       s[1]=s_ne[1];
       s[2]=s_ne[2];
       neighb_pt->object_pt()->get_x(s,x_large);

       //Get coordinates in small element
       s[0]=S_base(0,direction)+S_steplo(0,direction)+S_stephi(0,direction);
       s[1]=S_base(1,direction)+S_steplo(1,direction)+S_stephi(1,direction);
       s[2]=S_base(2,direction)+S_steplo(2,direction)+S_stephi(2,direction);
       el_pt->object_pt()->get_x(s,x_small);

       // Output
       if (neighbours_file.is_open())
        {
         neighbours_file << "#NORTH EAST: " << std::endl;
         neighbours_file << x_large[0]
                         << " " << x_large[1]
                         << " " << x_large[2]
                         << "  80 \n";
        }


       error=
        sqrt(pow(x_small[0]-x_large[0],2)+
             pow(x_small[1]-x_large[1],2)+
             pow(x_small[2]-x_large[2],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);
        }


       // Check error and doc mismatch if required
       if (std::fabs(error)>Max_neighbour_finding_tolerance) stop=true;
       if (stop)
        {
         if (!mismatch_file.is_open())
          {
           mismatch_file.open("mismatch.dat");
          }
         mismatch_file << "ZONE" << std::endl;
         mismatch_file << x_large[0] << " " << x_large[1] << " "
                       << x_large[2] <<"  2 \n";
         mismatch_file << x_small[0] << " " << x_small[1] << " "
                       << x_small[2] <<"  3 \n";
         mismatch_file.close();
         oomph_info << "Error for direction " << Direct_string[direction]
                    << " with element "<<i+1<< std::endl;
         pause("Error");
        }
      }

     // If neighbour does not exist: Insert blank zones into file
     // so that tecplot can find six neighbours for every element
     else
      {
       if (neighbours_file.is_open())
        {

         neighbours_file << "#---------------------------------" << std::endl;
         neighbours_file << "#no Neighbour in direct: "
                         << Direct_string[direction]
                         << std::endl;
         neighbours_file << "#---------------------------------" << std::endl;

         dynamic_cast<RefineableQElement<3>*>(
         el_pt->object_pt())->output_corners(neighbours_file,"WHITE");
         neighbours_file << "ZONE I=2 \n";
         neighbours_file << "-0.05 -0.05 -0.05  0 \n";
         neighbours_file << "-0.05 -0.05 -0.05  0 \n";
        }
      }

     if (neighbours_file.is_open())
      {
       neighbours_file << std::endl << std::endl << std::endl;
      }
    }
    } //End of case when element can be documented
  }
}


//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////



//=================================================================
/// Doc/check all true edge neighbours of octree (nodes) contained 
/// in the Vector forest_node_pt. Output into neighbours_file which can
/// be viewed from tecplot with OcTreeNeighbours.mcr
/// Neighbour info and errors are displayed on 
/// neighbours_txt_file.  Finally, compute the max. error between 
/// vertices when viewed from neighhbouring element. 
/// If the two filestreams are closed, output is suppressed. 
/// (Static function).
//=================================================================
void OcTree::doc_true_edge_neighbours(Vector<Tree*> forest_nodes_pt,
                                      std::ofstream& neighbours_file,
                                      std::ofstream& no_true_edge_file,
                                      std::ofstream& neighbours_txt_file,
                                      double& max_error)
{
 using namespace OcTreeNames;

 int diff_level;
 int edge=OMEGA;
 
 Vector<double> s(3);
 Vector<double> x(3);

 Vector<double> s_lo(3);
 Vector<double> s_hi(3);
 Vector<unsigned> translate_s(3);

 Vector<double> x_small(3);
 Vector<double> x_large(3);

 
 // Initialise error in vertex positions
 max_error=0.0;

 // Loop over all elements
 // ----------------------
 unsigned long num_nodes=forest_nodes_pt.size();

 for (unsigned long i=0;i<num_nodes;i++)
  {
   
   // Doc the element itself
   OcTree* el_pt = dynamic_cast<OcTree*>(forest_nodes_pt[i]);

   //If the object is incomplete omit
   if(el_pt->object_pt()->nodes_built())
    {
     // Print it
     if (neighbours_file.is_open())
      {
       neighbours_file << "#---------------------------------" << std::endl;
       neighbours_file << "#The element itself: " << i << std::endl;
       neighbours_file << "#---------------------------------" << std::endl;
       dynamic_cast<RefineableQElement<3>*>(
        el_pt->object_pt())->output_corners(neighbours_file,"BLACK");
      }

   // Loop over directions to find edge neighbours
   // --------------------------------------------
   for (int direction=LB;direction<=UF;direction++)
    {
     // Initialise difference in levels 
     diff_level=0;

     // For now simply doc the zero-th edge neighbour (if any)
     unsigned i_root_edge_neighbour=0;
     unsigned nroot_edge_neighbour=0; 

     // Find greater-or-equal-sized edge neighbour...
     OcTree* neighb_pt=el_pt->gteq_true_edge_neighbour(direction,
                                                       i_root_edge_neighbour,
                                                       nroot_edge_neighbour, 
                                                       translate_s,
                                                       s_lo, s_hi,
                                                       edge,diff_level);
     
      // If neighbour exists and nodes are created
      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
                            << ". Inside the neighbour the edge is "
                            << Direct_string[edge] <<  std::endl;
        }

       // Plot neighbour in the appropriate colour
       if (neighbours_file.is_open())
        {
         neighbours_file << "#---------------------------------" << std::endl;
         neighbours_file << "#Neighbour element: " << Direct_string[direction]
                         << std::endl;
         neighbours_file << "#---------------------------------" << std::endl;
         dynamic_cast<RefineableQElement<3>*>(
          neighb_pt->object_pt())->
          output_corners(neighbours_file,Colour[direction]);
        }

       // Check that local coordinates in the larger element
       // lead to the same spatial point as the node vertices
       // in the current element
       if (neighbours_file.is_open())
        {
         neighbours_file << "ZONE I=2  C=" << Colour[direction] << std::endl;
        }

       // "Low" vertex in the edge
       //-------------------------

       // Get coordinates in large (neighbour) element
       s[0]=s_lo[0];
       s[1]=s_lo[1];
       s[2]=s_lo[2];
       neighb_pt->object_pt()->get_x(s,x_large);

       // Get coordinates in small element
       Vector<double> s(3);
       s[0]=S_base_edge(0,direction);
       s[1]=S_base_edge(1,direction);
       s[2]=S_base_edge(2,direction);
       el_pt->object_pt()->get_x(s,x_small);

       double error=
        sqrt(pow(x_small[0]-x_large[0],2)+
             pow(x_small[1]-x_large[1],2)+
             pow(x_small[2]-x_large[2],2));

       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);
        }


       // Check error and doc mismatch if required
       bool stop=false;
       std::ofstream mismatch_file;
       if (std::fabs(error)>Max_neighbour_finding_tolerance)
        {
         stop=true;
         mismatch_file.open("mismatch.dat");
         mismatch_file << "ZONE" << std::endl;
         mismatch_file<< x_large[0]<<" "<<x_large[1]<<" "<<x_large[2]<<" 2 \n";
         mismatch_file<< x_small[0]<<" "<<x_small[1]<<" "<<x_small[2]<<" 3 \n";
        }

       if (neighbours_file.is_open())
        {
         neighbours_file << "#LOW VERTEX: " << std::endl;
         neighbours_file << x_large[0] << " "
                         << x_large[1] << " "
                         << x_large[2] <<"  40 \n";
        }


       // "Hi" vertex in the edge
       //------------------------

       // Get coordinates in large (neighbour) element
       s[0]=s_hi[0];
       s[1]=s_hi[1];
       s[2]=s_hi[2];
       neighb_pt->object_pt()->get_x(s,x_large);

       //Get coordinates in small element
       s[0]=S_base_edge(0,direction)+S_step_edge(0,direction);
       s[1]=S_base_edge(1,direction)+S_step_edge(1,direction);
       s[2]=S_base_edge(2,direction)+S_step_edge(2,direction);
       el_pt->object_pt()->get_x(s,x_small);

       // Output
       if (neighbours_file.is_open())
        {
         neighbours_file << "#HI VERTEX: " << std::endl;
         neighbours_file << x_large[0]
                         << " " << x_large[1]
                         << " " << x_large[2]
                         << "  80 \n";
        }


       error=
        sqrt(pow(x_small[0]-x_large[0],2)+
             pow(x_small[1]-x_large[1],2)+
             pow(x_small[2]-x_large[2],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);
        }


       // Check error and doc mismatch if required
       if (std::fabs(error)>Max_neighbour_finding_tolerance) stop=true;
       if (stop)
        {
         if (!mismatch_file.is_open())
          {
           mismatch_file.open("mismatch.dat");
          }
         mismatch_file << "ZONE" << std::endl;
         mismatch_file << x_large[0] << " " << x_large[1] << " "
                       << x_large[2] <<"  2 \n";
         mismatch_file << x_small[0] << " " << x_small[1] << " "
                       << x_small[2] <<"  3 \n";
         mismatch_file.close();
         oomph_info << "Error for direction " << Direct_string[direction]
                    << " with element "<<i+1<< std::endl;
         pause("Error");
        }
      }

     // If neighbour does not exist: Insert blank zones into file
     // so that tecplot can find twelve neighbours for every element
     else
      {
       if (neighbours_file.is_open())
        {

         neighbours_file << "#---------------------------------" << std::endl;
         neighbours_file << "#no Neighbour in direct: "
                         << Direct_string[direction]
                         << std::endl;
         neighbours_file << "#---------------------------------" << std::endl;

         dynamic_cast<RefineableQElement<3>*>(
         el_pt->object_pt())->output_corners(neighbours_file,"WHITE");
         neighbours_file << "ZONE I=2 \n";
         neighbours_file << "-0.05 -0.05 -0.05  0 \n";
         neighbours_file << "-0.05 -0.05 -0.05  0 \n";



         // Doc edge for which no neighbour exists but only
         // for the smallest elements
         if (el_pt->is_leaf())
          {
           // Get start coordinates in small element
           Vector<double> s(3),x_start(3),x_end(3);
           s[0]=S_base_edge(0,direction);
           s[1]=S_base_edge(1,direction);
           s[2]=S_base_edge(2,direction);
           el_pt->object_pt()->get_x(s,x_start);
           
           
           //Get coordinates in small element
           s[0]=S_base_edge(0,direction)+S_step_edge(0,direction);
           s[1]=S_base_edge(1,direction)+S_step_edge(1,direction);
           s[2]=S_base_edge(2,direction)+S_step_edge(2,direction);
           el_pt->object_pt()->get_x(s,x_end);
           
           no_true_edge_file << "ZONE I=2" << std::endl;
           no_true_edge_file << x_start[0] << " "
                             << x_start[1] << " "
                             << x_start[2] << " "
                             << std::endl;
           no_true_edge_file << x_end[0] << " "
                             << x_end[1] << " "
                             << x_end[2] << " "
                             << std::endl;
           
          }
        }
       
       if (neighbours_file.is_open())
        {
         neighbours_file << std::endl << std::endl << std::endl;
        }
      }
    } //End of case when element can be documented
   
    }
  }
  
}



//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////



//================================================================
/// Constructor for OcTreeForest:
///
/// Pass:
///  - trees_pt[], the Vector of pointers to the constituent trees
///    (OcTreeRoot objects)
///
//=================================================================
OcTreeForest::OcTreeForest(Vector<TreeRoot* >& trees_pt) :
 TreeForest(trees_pt)
{
#ifdef LEAK_CHECK
 LeakCheckNames::OcTreeForest_build+=1;
#endif

 // Don't setup neighbours etc. if forest is empty
 if (trees_pt.size()==0)
  {
   return;
  }

 using namespace OcTreeNames;
 
 //Construct the neighbour and rotation scheme, note that all neighbour 
 //pointers must be set before the constructor is called
 
//  MemoryUsage::doc_memory_usage(
//   "before find_neighbours in octree forest constr");

 // setup the neighbour scheme
 find_neighbours();

//  MemoryUsage::doc_memory_usage(
//   "after find_neighbours in octree forest constr");
 
 // setup the rotation scheme
 construct_up_right_equivalents();

//  MemoryUsage::doc_memory_usage(
//   "after construct_up_right_equivalents in octree forest constr");

}

//================================================================
/// setup the neighbour scheme : tells to each element in the 
/// forest which (if any) element is its {R,L,U,D,B,F} face neighbour
/// and which is its {LB,RB,...,UF} edge neighbour.
//================================================================
void OcTreeForest::find_neighbours()
{
 using namespace OcTreeNames ;
 unsigned numtrees=ntree();
 int 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",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }


 // Number of vertex nodes: 8
 unsigned n_vertex_node=8;

 // 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<BrickElementBase*>(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++)
  {

   // Cast to OcTreeRoot
   OcTreeRoot* octree_root_i_pt=dynamic_cast<OcTreeRoot*>(Trees_pt[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);

     // So far, we haven't identified the candidate element as a face
     // neighbour of element/tree i
     bool is_face_neighbour=false;

     // is it the Up neighbour ?
     bool is_Up_neighbour=
      ((Trees_pt[j]->object_pt()->get_node_number(
         Trees_pt[i]->object_pt()->node_pt(n*(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 Down neighbour ?
     bool is_Down_neighbour=
      ((Trees_pt[j]->object_pt()->get_node_number(
         Trees_pt[i]->object_pt()->node_pt(n*n*(n-1)))!=-1) &&
       (Trees_pt[j]->object_pt()->get_node_number(
        Trees_pt[i]->object_pt()->node_pt(n-1))!=-1));
      
     // is it the Right neighbour ?
     bool is_Right_neighbour=
      ((Trees_pt[j]->object_pt()->get_node_number(
         Trees_pt[i]->object_pt()->node_pt(n*n*n-1))!=-1) &&
       (Trees_pt[j]->object_pt()->get_node_number(
        Trees_pt[i]->object_pt()->node_pt(n-1))!=-1));
       
     // is it the Left neighbour ?
     bool is_Left_neighbour=
      ((Trees_pt[j]->object_pt()->get_node_number(
         Trees_pt[i]->object_pt()->node_pt(n*(n*n-1)))!=-1) &&
       (Trees_pt[j]->object_pt()->get_node_number(
        Trees_pt[i]->object_pt()->node_pt(0))!=-1));
 
     // is it the Back neighbour ?
     bool is_Back_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(0))!=-1));

     // is it the Front neighbour ?
     bool is_Front_neighbour=
      ((Trees_pt[j]->object_pt()->get_node_number(
         Trees_pt[i]->object_pt()->node_pt(n*n*n-1))!=-1) &&
       (Trees_pt[j]->object_pt()->get_node_number(
        Trees_pt[i]->object_pt()->node_pt(n*n*(n-1)))!=-1));
      

     if(is_Down_neighbour)
      {
       is_face_neighbour=true;
       Trees_pt[i]->neighbour_pt(D)=Trees_pt[j];
      }
     if(is_Up_neighbour)
      {
       is_face_neighbour=true; 
       Trees_pt[i]->neighbour_pt(U)=Trees_pt[j];
      }
     if(is_Right_neighbour)
      {
       is_face_neighbour=true; 
       Trees_pt[i]->neighbour_pt(R)=Trees_pt[j];
      }
     if(is_Left_neighbour)
      {
       is_face_neighbour=true; 
       Trees_pt[i]->neighbour_pt(L)=Trees_pt[j];
      }
     if(is_Back_neighbour)
      {
       is_face_neighbour=true; 
       Trees_pt[i]->neighbour_pt(B)=Trees_pt[j];
      }
     if(is_Front_neighbour)
      {
       is_face_neighbour=true; 
       Trees_pt[i]->neighbour_pt(F)=Trees_pt[j];
      }


     // If it's not a face neighbour, it may still be an edge
     // neighbour. We check this by checking if the
     // vertex nodes coincide. Note: This test would also
     // evaluate to true for face neighbours but we've already
     // determined that the element is not a face neighbour!
     if (!is_face_neighbour)
      {
         
       // is it the left back neighbour ?
       bool is_left_back_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));

       // is it the right back neighbour ?
       bool is_right_back_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 down back neighbour ?
       bool is_down_back_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(0))!=-1));

       // is it the up back neighbour ?
       bool is_up_back_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 left down neighbour ?
       bool is_left_down_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt(n*n*(n-1)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt(0))!=-1));


       // is it the right down neighbour ?
       bool is_right_down_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt((n*n+1)*(n-1)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt(n-1))!=-1));


       // is it the left up neighbour ?
       bool is_left_up_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt((n*n*n-n)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt(n*(n-1)))!=-1));


       // is it the right up neighbour ?
       bool is_right_up_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt((n*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 left front neighbour ?
       bool is_left_front_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt((n*n*n-n)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt(n*n*(n-1)))!=-1));


       // is it the right front neighbour ?
       bool is_right_front_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt((n*n*n-1)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt((n*n+1)*(n-1)))!=-1));
                  

       // is it the down front neighbour ?
       bool is_down_front_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt(n*n*(n-1)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt((n*n+1)*(n-1)))!=-1));


       // is it the up front neighbour ?
       bool is_up_front_neighbour=
        ((Trees_pt[j]->object_pt()->get_node_number(
           Trees_pt[i]->object_pt()->node_pt((n*n*n-n)))!=-1) &&
         (Trees_pt[j]->object_pt()->get_node_number(
          Trees_pt[i]->object_pt()->node_pt((n*n*n-1)))!=-1));


                  
       // Add to storage scheme for edge neighbours (only!)
       
       if(is_left_back_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(LB)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],LB);
        }
       if(is_right_back_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(RB)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],RB);
        }
       if(is_down_back_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(DB)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],DB);
        }
       if(is_up_back_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(UB)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],UB);
        }
         
         
       if(is_left_down_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(LD)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],LD);
        }
       if(is_right_down_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(RD)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],RD);
        }
       if(is_left_up_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(LU)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],LU);
        }
       if(is_right_up_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(RU)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],RU);
        }
         
         
       if(is_left_front_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(LF)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],LF);
        }
       if(is_right_front_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(RF)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],RF);
        }
       if(is_down_front_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(DF)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],DF);
        }
       if(is_up_front_neighbour)
        {
         //Trees_pt[i]->neighbour_pt(UF)=Trees_pt[j];
         octree_root_i_pt->add_edge_neighbour_pt(Trees_pt[j],UF);
        }

      }
    }
  }

}
   
 
 
//================================================================
/// Construct the rotation scheme for the octree forest.
/// Note that all pointers to neighbours must have been allocated
/// for this to work.
//================================================================
void OcTreeForest::construct_up_right_equivalents()
{
 using namespace OcTreeNames;   

 // Number of trees in forest
 unsigned numtrees = ntree();

 // nnode1d from first element (if it exists!)
 int nnode1d=0;
 if(numtrees>0) nnode1d=Trees_pt[0]->object_pt()->nnode_1d();

 OcTreeRoot* neigh_pt=0;
 


 //Loop over all the trees
 //------------------------
 for(unsigned i=0;i<numtrees;i++)
  {

   //Find the pointer to the Up neighbour
   //------------------------------------
   neigh_pt = oc_face_neigh_pt(i,U);

   //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(octree_pt(i));

     //If up neighbour has a pointer to this tree
     if(direction!=Tree::OMEGA)
      {
       //The direction to the element in the neighbour 
       //must be equivalent to the down direction in this element
       //Hence, the up equivalent is the reflection of that direction
       octree_pt(i)->set_up_equivalent(neigh_pt,
                                       OcTree::Reflect_face[direction]);

       //The right equivalent is the direction to the equivalent of
       //the right face in the present element, which will
       //be connected to the UR edge of the present element,
       //but will not be the face adjacent to the U boundary
       //i.e. the "direction" face.
       //We find the local node numbers, in the neighbour, of
       //the UR edge of the present element 
       int nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d*nnode1d-1));
       int nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d-1));
                                                       
       //Now get the other face connected to that edge in the
       //neighbour. It is the right equivalent
       octree_pt(i)->set_right_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
      }
     else     
      //If U neighbour does not have pointer to this tree, die
      {
       std::ostringstream error_stream;
       error_stream 
        << "Tree " << i
        << "'s Up 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 Right neighbour
   //---------------------------------------
   neigh_pt = oc_face_neigh_pt(i,R);

   //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(octree_pt(i));

     //If the neighbour has a pointer to this tree
     if(direction != Tree::OMEGA)
      {

       //The direction to the element in the neighbour
       //must be equivalent to the left direction in this element
       //Hence, the right equivalent is the reflection of that direction
       octree_pt(i)->set_right_equivalent(neigh_pt,
                                          OcTree::Reflect_face[direction]);
       
       //The up equivalent will be the direction to the equivalent of
       //the up face in the neighbour, which will be connected to the
       //UR edge of the present element, but will not be the face 
       //adjacent to the R boundary, i.e. the "direction" face
       //We find the local node numbers, in the neighbour, of the
       //UR edge of the present element
       int nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d*nnode1d-1));
       int nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d-1));
   
       //Now get the other face connected to that edge in the 
       //neighbour. It is the up equivalent
       octree_pt(i)->set_up_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));

      }
     else
      {
       //If R neighbour does not have pointer to this tree, die
       std::ostringstream error_stream;
       error_stream 
        << "Tree " << i
        << "'s Right 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 Down neighbour
   //--------------------------------------
   neigh_pt = oc_face_neigh_pt(i,D);

   //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(octree_pt(i));

     //If the neighbour has a pointer to this element
     if(direction!=Tree::OMEGA)
      {
       //The direction to the element in the neighbour must be
       //equivalent to the up direction
       octree_pt(i)->set_up_equivalent(neigh_pt,direction);

       //The right equivalent is the direction to the equivalent of
       //the right face in the present element, which will be 
       //connected to the DR edge of the present element, but will
       //not be the face adjacent to the D boundary,
       //i.e. the "direction" face.
       //We find the local node numbers, in the neighbour, of
       //the RD edge of the present element
       int nod1=neigh_pt->object_pt()
        ->get_node_number(octree_pt(i)->object_pt()
                          ->node_pt(nnode1d-1));
       int nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt((nnode1d*nnode1d+1)*(nnode1d-1)));

       //Now get the other face connected to that edge in the neighbour.
       //It is the right equivalent
       octree_pt(i)->set_right_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
      }    
     else
      {
       //If D neighbour does not have pointer to this tree, die
       std::ostringstream error_stream;
       error_stream 
        << "Tree " << i
        << "'s Down 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 Left neighbour
   //--------------------------------------
   neigh_pt = oc_face_neigh_pt(i,L);

   //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(octree_pt(i));

     //If the neighbour has a pointer to the present element
     if(direction!=Tree::OMEGA)
      {
       //The direction to the element in the neighbour is
       //must be equivalent to the right direction
       octree_pt(i)->set_right_equivalent(neigh_pt,direction);
       
       //The up equivalent is the direction to the equivalent of the
       //up face in the present element, which will be connected to
       //the UL edge of the present element, but will not
       //be the face adjacent to the L boundary, i.e. the "direction"
       //face. 
       //We find the local node numbers, in the neighbour, of the UL
       //edge
       int nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt((nnode1d-1)*nnode1d));
       int nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt((nnode1d*nnode1d-1)*nnode1d));
  
       //Now get the other face connected to that edge in the
       //neighbour.It is the up equivalent
       octree_pt(i)->set_up_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
       
      }
     else
      {
       //If L neighbour does not have pointer to this tree, die
       std::ostringstream error_stream;
       error_stream 
        << "Tree " << i
        << "'s Left 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 Front neighbour
   //---------------------------------------
   neigh_pt = oc_face_neigh_pt(i,F);

   //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(octree_pt(i));

     //If the neighbour has a pointer to the present element
     if(direction!=Tree::OMEGA)
      {
       //The direction to the up face will be given by the
       //other face connected to the UF edge in the neighbour
       //i.e. the face that is not given by direction
       //We obtain the local node numbers, in the neighbour,
       //of the UF edge in the present element
       int nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d*nnode1d-1));
       int nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt((nnode1d*nnode1d-1)*nnode1d));

       //We now get the other face connected to that edge
       //It is the up equivalent.
       octree_pt(i)->set_up_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
       
       //The direction to the right face will be given by the
       //other face connected to the RF edge in the neighbour
       //i.e. the face that is not given by the direction
       //We get the local node numbers, in the neighbour,
       //of the RF edge in the present element
       nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d*nnode1d-1));
       nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt((nnode1d*nnode1d+1)*(nnode1d-1)));

       //We now get the other face connected to that edge
       //It is the right equivalent
       octree_pt(i)->set_right_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
      }
     else
      {
       //If F neighbour does not have pointer to this tree, die
       std::ostringstream error_stream;
       error_stream 
        << "Tree " << i
        << "'s Front 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 Back neighbour
   //--------------------------------------
   neigh_pt = oc_face_neigh_pt(i,B);

   //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(octree_pt(i));

     //If the neighbour has a pointer to the present element
     if(direction!=Tree::OMEGA)
      {
       //The direction to the up face will be given by the
       //other face connected to the UB edge of the present element
       //i.e. not the "direction" face
       //We obtain the local node numbers, in the neighbour,
       //of the UB edge in the present element
       int nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt((nnode1d-1)*nnode1d));
       int nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d-1));

       //We now get the other face connected to that edge
       //It is the up equivalent
       octree_pt(i)->set_up_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
       
       
       //The direction to the right face will be given by
       //the other face connected to the RB edge of the present
       //element, i.e. not the direction face.
       //We obtain local node numbers, in the neighbour,
       //of the RB edge in the present element
       nod1=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d*nnode1d-1));
       nod2=neigh_pt->object_pt()->
        get_node_number(octree_pt(i)->object_pt()
                        ->node_pt(nnode1d-1));

       //We now get the other face connected to that edge
       //It is the right equivalent
       octree_pt(i)->set_right_equivalent(
        neigh_pt,
        OcTree::get_the_other_face(nod1,nod2,nnode1d,direction));
      }
     else
      {
       //If B neighbour does not have pointer to this tree, die
       std::ostringstream error_stream;
       error_stream 
        << "Tree " << i
        << "'s Back neighbour has no neighbour pointer to Tree " << i
        << std::endl;

       throw OomphLibError(error_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }


   // EDGE NEIGHBOURS
   //----------------

   // Loop over all edges:
   for (int edge=LB;edge<=UF;edge++)
    {
     // Get vector to pointers of edge neighbours
     Vector<TreeRoot*> edge_neigh_pt=oc_edge_neigh_pt(i,edge);
     
     // Loop over all edge neighbours
     unsigned n_neigh=edge_neigh_pt.size();
     for (unsigned e=0;e<n_neigh;e++)
      {
       //If there is a neighbour
       if(edge_neigh_pt[e]!=0)
        {         
         // Here are the two vertices at the two ends of the edge
         // in the present element
         int vertex      =OcTree::Vertex_at_end_of_edge[edge][0];       
         int other_vertex=OcTree::Vertex_at_end_of_edge[edge][1];       
         
         // Get node numbers of the two vertices on edge
         unsigned nod1=OcTree::vertex_to_node_number(vertex,nnode1d);       
         unsigned nod2=OcTree::vertex_to_node_number(other_vertex,nnode1d);
         
         
         // Here are the local node numbers (in the neighbouring element)
         // of the start/end nodes of present element's edge:
         unsigned neighb_nod1=edge_neigh_pt[e]->object_pt()->
          get_node_number(octree_pt(i)->object_pt()
                          ->node_pt(nod1));
         
         unsigned neighb_nod2=edge_neigh_pt[e]->object_pt()->
          get_node_number(octree_pt(i)->object_pt()
                          ->node_pt(nod2));
         
         // Convert to vertices
         int neighb_vertex=
          OcTree::node_number_to_vertex(neighb_nod1,nnode1d);
         int neighb_other_vertex=
          OcTree::node_number_to_vertex(neighb_nod2,nnode1d);
                  
         //Up equivalent is stored first in the pair that's returned
         //by the lookup table
         octree_pt(i)->set_up_equivalent(
          edge_neigh_pt[e],
          OcTree::Up_and_right_equivalent_for_pairs_of_vertices[
           std::make_pair(std::make_pair(neighb_vertex,vertex),
                          std::make_pair(neighb_other_vertex,other_vertex))].
          first);

         //Right equivalent is stored second in the pair that's returned
         //by the lookup table
         octree_pt(i)->set_right_equivalent(edge_neigh_pt[e],
          OcTree::Up_and_right_equivalent_for_pairs_of_vertices[
           std::make_pair(std::make_pair(neighb_vertex,vertex),
                          std::make_pair(neighb_other_vertex,other_vertex))].
                                            second);
         
        }
      }
    }
  } //end of loop over all trees.
}




//================================================================
/// Self test: Check neighbour finding routine
//=================================================================
unsigned OcTreeForest::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 vertices and their opposite points
 std::ofstream neighbours_file;
 std::ofstream no_true_edge_file;
 std::ofstream neighbours_txt_file;
 
 double max_error_face=0.0;
 OcTree::doc_face_neighbours(all_nodes_pt,neighbours_file,
                             neighbours_txt_file, max_error_face);
 
 double max_error_edge=0.0;
 OcTree::doc_true_edge_neighbours(all_nodes_pt,neighbours_file,
                                  no_true_edge_file,
                                  neighbours_txt_file, max_error_edge);
 
 bool failed=false;
 if (max_error_face>OcTree::max_neighbour_finding_tolerance())
  {
   oomph_info 
    << "\n\n Failed self_test() for OcTreeForest because of faces: Max. error "
    << max_error_face << std::endl<< std::endl;
   failed=true;
  }

 if (max_error_edge>OcTree::max_neighbour_finding_tolerance())
  {
   oomph_info 
    << "\n\n Failed self_test() for OcTreeForest because of edges: Max. error "
    << max_error_edge << std::endl<< std::endl;
   failed=true;
  }

 double max_error=max_error_face;
 if (max_error_edge>max_error) max_error=max_error_edge;

 if (failed)
  {
   return 1;
  }
 else
  {
   oomph_info << "\nPassed self_test() for OcTreeForest: Max. error "
              << max_error << std::endl;
   return 0;
  }
 
}




//================================================================
/// Document and check all the neighbours of all the nodes
/// in the forest. DocInfo object specifies the output directory
/// and file numbers for the various files. If 
/// \c doc_info.is_doc_enabled()=false
/// no output is created.
//================================================================
void OcTreeForest::check_all_neighbours(DocInfo &doc_info)
{
 // Create vector of all elements in the tree
 Vector<Tree*> all_tree_nodes_pt;
 this->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);

 // Face neighbours
 //----------------
 {
  //Create storage for the files
  std::ofstream neigh_file;
  std::ofstream neigh_txt_file;
  //If we are documenting, then do so
  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 static member of OcTree function
  double max_error=0.0;
  OcTree::doc_face_neighbours(all_tree_nodes_pt,
                              neigh_file,neigh_txt_file,max_error);
  if (max_error>Tree::max_neighbour_finding_tolerance())
   {
    std::ostringstream error_stream;
    error_stream << "\nMax. error in octree 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 << "\nMax. error in octree neighbour finding: " 
               << max_error << " is OK" << std::endl;
    oomph_info << "i.e. less than OcTree::max_neighbour_finding_tolerance()=" 
               << OcTree::max_neighbour_finding_tolerance() << std::endl;
   }
  
  //Close the documentation files if they were opened
  if(doc_info.is_doc_enabled())
   {
    neigh_file.close();
    neigh_txt_file.close();
   }
 }

 // Edge neighbours
 //----------------
 {
  //Create storage for the files
  std::ofstream neigh_file;
  std::ofstream no_true_edge_file;
  std::ofstream neigh_txt_file;
  //If we are documenting, then do so
  if (doc_info.is_doc_enabled())
   {
    std::ostringstream fullname;
    fullname << doc_info.directory() << "/edge_neighbours"
             << doc_info.number() << ".dat";
    neigh_file.open(fullname.str().c_str());
    fullname.str("");
    fullname << doc_info.directory() << "/no_true_edge"
             << doc_info.number() << ".dat";
    no_true_edge_file.open(fullname.str().c_str());
    fullname.str("");
    fullname << doc_info.directory() << "/edge_neighbours"
             << doc_info.number() << ".txt";
    neigh_txt_file.open(fullname.str().c_str());  
   }
  
  //Call the static member of OcTree function
  double max_error=0.0;
  OcTree::doc_true_edge_neighbours(all_tree_nodes_pt,
                                   neigh_file,no_true_edge_file,
                                   neigh_txt_file,max_error);
  if (max_error>Tree::max_neighbour_finding_tolerance())
   {
    std::ostringstream error_stream;
    error_stream << "Max. error in octree edge 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();
      no_true_edge_file.close();
      neigh_txt_file.close();
     }

    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  else
   {
    oomph_info << "Max. error in octree edge neighbour finding: " 
               << max_error << " is OK" << std::endl;
    oomph_info << "i.e. less than OcTree::max_neighbour_finding_tolerance()=" 
               << OcTree::max_neighbour_finding_tolerance() << std::endl;
   }
  
  //Close the documentation files if they were opened
  if(doc_info.is_doc_enabled())
   {
    neigh_file.close();
    no_true_edge_file.close();
    neigh_txt_file.close();
   }
 }


}



//================================================================
/// Open output files that will stored any hanging nodes that are
/// created in the mesh refinement process.
///===============================================================
void OcTreeForest::open_hanging_node_files(DocInfo &doc_info,
                                           Vector<std::ofstream*> 
                                           &output_stream)
{
 //In 3D, there will be six output files
 for(unsigned i=0;i<6;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_u"
            << doc_info.number() << ".dat";
   output_stream[0]->open(fullname.str().c_str());
   fullname.str("");
   fullname << doc_info.directory() << "/hang_nodes_d"
            << doc_info.number() << ".dat";
   output_stream[1]->open(fullname.str().c_str());
   fullname.str("");
   fullname << doc_info.directory() << "/hang_nodes_l"
            << doc_info.number() << ".dat";
   output_stream[2]->open(fullname.str().c_str());
   fullname.str("");
   fullname << doc_info.directory() << "/hang_nodes_r"
            << doc_info.number() << ".dat";
   output_stream[3]->open(fullname.str().c_str());
   fullname.str("");
   fullname << doc_info.directory() << "/hang_nodes_b"
            << doc_info.number() << ".dat";
   output_stream[4]->open(fullname.str().c_str());
   fullname.str("");
   fullname << doc_info.directory() << "/hang_nodes_f"
            << doc_info.number() << ".dat";
   output_stream[5]->open(fullname.str().c_str());
  }
}



}