octree.h

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//====================================================================
#ifndef OOMPH_OCTREE_HEADER
#define OOMPH_OCTREE_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

//OOMPH-LIB headers
#include "tree.h"
#include "matrices.h"

namespace oomph
{



//====================================================================
// Namespace for OcTree directions
//====================================================================
namespace OcTreeNames
{

 /// \short Directions. OMEGA is used if a direction is undefined 
 /// in a certain context
 enum{ 
       LDB,RDB,LUB,RUB,     // back octants/vertices
       LDF,RDF,LUF,RUF,     // front octants/vertices
       LB,RB,DB,UB,         // back edges
       LD,RD,LU,RU,         // side edges
       LF,RF,DF,UF,         // front edges
       L,R,D,U,B,F,          // faces
       OMEGA=26};
}


// Forward definitions
class OcTreeRoot;

//================================================================
/// OcTree class: Recursively defined, generalised octree.
///
/// An OcTree has:
/// - a pointer to the object (of type  RefineableQElement<3>) it represents 
/// - Vector of pointers to its eight (LDB,RDB,...,RUF) sons (which are 
///   octrees themselves). If the Vector of pointers to the sons has zero
///   length, the OcTree is a "leaf node" in the overall octree.
/// - a pointer to its father. If this pointer is NULL, the OcTree is the
///   the root node of the overall octree.
/// This data is stored in the Tree base class. 
/// 
/// The tree can also be part of a forest. If that is the case, the root
/// will have pointers to the roots of neighbouring octrees.
/// 
/// The objects contained in the octree are assumed to be
/// (topologically) cubic elements whose geometry is
/// parametrised by local coordinates \f$ {\bf s} \in [-1,1]^3 \f$.
///
/// The tree can be traversed while actions are being performed
/// at all of its "nodes" or only at the leaf "nodes".
///
/// Finally, the leaf "nodes" can be split depending on criteria
/// defined by the object.
///
/// Note that OcTrees are only generated by splitting existing
/// OcTrees. Therefore, the constructors are protected. The
/// only OcTree that "Joe User" can create is
/// the (derived) class OcTreeRoot.
//=================================================================
class OcTree : public virtual Tree
{ 

 public:

  /// \short Destructor. Note: Deleting a octree also deletes the 
  /// objects associated with all non-leaf nodes!
 virtual ~OcTree() {}


 /// Broken copy constructor
 OcTree(const OcTree& dummy) 
  { 
   BrokenCopy::broken_copy("OcTree");
  } 
 
 /// Broken assignment operator
 void operator=(const OcTree&) 
  {
   BrokenCopy::broken_assign("OcTree");
  }


 /// \short Overload the function construct_son to ensure that the son
 /// is a specific OcTree and not a general Tree.
 Tree* construct_son(RefineableElement* const &object_pt,
                     Tree* const &father_pt, const int &son_type)
  {
   OcTree* temp_oc_pt = new OcTree(object_pt,father_pt,son_type);
   return temp_oc_pt;
  }

 
/// \short 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* gteq_face_neighbour(const int& direction,
                             Vector<unsigned> &translate_s,
                             Vector<double>& s_sw,
                             Vector<double>& s_ne, 
                             int& face,
                             int& diff_level) const;



/// \short 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* 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;
 

 /// \short Self-test: Check all neighbours. Return success (0) 
 /// if the max. distance between corresponding points in the
 /// neighbours is less than the tolerance specified in the
 /// static value Tree::Max_neighbour_finding_tolerance.
 unsigned self_test();

 /// \short Setup the static data, rotation and reflection schemes, etc
 static void setup_static_data();


 /// \short 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 void doc_face_neighbours(Vector<Tree*> forest_nodes_pt, 
                                 std::ofstream& neighbours_file, 
                                 std::ofstream& neighbours_txt_file,
                                 double& max_error);


 /// \short 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 void 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);
 
 /// \short If an edge is bordered by the nodes whose local numbers 
 /// are n1 and n2 in an element with nnode1d nodes along each coordinate
 /// direction, then this edge is shared by two faces. This function 
 /// takes one of these faces as the argument \c face and returns the 
 /// other one. (\c face is a direction in the set U,D,F,B,L,R).
 static int get_the_other_face(const unsigned& n1, const unsigned& n2, 
                               const unsigned& nnode1d, const int& face); 

 /// \short Return the local node number of given vertex
 /// [LDB,RDB,...] in an element with nnode1d nodes in each 
 /// coordinate direction
 static unsigned vertex_to_node_number(const int& vertex, 
                                       const unsigned& nnode1d);


  /// \short Return the vertex  [LDB,RDB,...] of local (vertex) node n
 /// in an element with nnode1d nodes in each coordinate direction.
 static int node_number_to_vertex(const unsigned& n, 
                                  const unsigned& nnode1d);


 /// \short If U[p] becomes new_up and R[ight] becomes new_right then the 
 /// direction vector \c dir becomes rotate(new_up, new_right, dir)
 static Vector<int> rotate(const int& new_up, const int& new_right,
                           const Vector<int>& dir);


 /// \short If U[p] becomes new_up and R[ight] becomes new_right 
 /// then the direction \c dir becomes \c rotate(new_up, new_right, dir)
 static int rotate(const int& new_up, const int& new_right,
                   const int& dir);



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

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

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

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

 /// \short \c Vector of vectors containing the two vertices for each edge,
 /// e.g. \c Vertex_at_end_of_edge[LU][0]=LUB and
 /// \c Vertex_at_end_of_edge[LU][1]=LUF.
 static Vector<Vector<int> > Vertex_at_end_of_edge;

 /// \short Each vector representing a direction can be translated into
 /// a direction, either a son type (vertex), a face or an edge. 
 /// E.g. : Vector_to_direction[(1,-1,1)]=RDF, Vector_to_direction[(0,1,0)]=U  
 static std::map<Vector<int>, int > Vector_to_direction;
 
 /// \short For each direction, i.e. a son_type (vertex), a face or an 
 /// edge, this defines a vector that indicates this direction.
 /// E.g : Direction_to_vector[RDB]=(1,-1,-1), Direction_to_vector[U]=(0,1,0)
 static Vector<Vector<int> > 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. 
 static std::map<std::pair<std::pair<int,int>,std::pair<int,int> >, 
  std::pair<int,int> > Up_and_right_equivalent_for_pairs_of_vertices;


  protected:

 /// Default constructor (empty and broken)
 OcTree() 
  {
   throw OomphLibError(
    "Don't call empty constructor for OcTree!",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Constructor for empty (root) tree: 
 /// no father, no sons; just pass a pointer to its object
 /// (a RefineableQElement<3>). This is
 /// protected because OcTrees can only be created internally,
 /// during the split operation. Only OcTreeRoots can be
 /// created externally. 
 OcTree(RefineableElement* const &object_pt) : Tree(object_pt) {}


 /// \short Constructor for tree that has a father: Pass it the pointer 
 /// to its object, the pointer to its father and tell it what type 
 /// of son (LDB,RDB,...) it is.
 /// Protected because OcTrees can only be created internally,
 /// during the split operation.  Only OcTreeRoots can be
 /// created externally. 
 OcTree(RefineableElement* const &object_pt, 
        Tree* const &father_pt, const int& son_type) :
  Tree(object_pt,father_pt,son_type) {}

 
 /// Bool indicating that static member data has been setup
 static bool Static_data_has_been_setup;


  private:

 /// \short 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  the maximum level to which 
 /// search is performed because in the presence of OcTree forests,
 /// the search isn't purely recursive.
 /// 
 /// 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.
 /// - 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* gteq_face_neighbour(const int& direction, 
                             double& s_difflo, 
                             double& s_diffhi,
                             int& diff_level ,
                             int max_level, 
                             OcTreeRoot* orig_root_pt) const;


 /// \short 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  the maximum level to which 
 /// search is performed because in the presence of OcTree forests,
 /// the search isn't purely recursive.
 /// 
 /// 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 the edge's "low" 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.
 /// .
 /// \b Note: some of the auxiliary information may be incorrect if
 /// the neighbour is not a true edge neighbour.  We don't care because
 /// we're not dealing with those!
 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;
 

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




 /// \short This constructs the rotation matrix of the rotation around the 
 /// axis \c axis with an angle of \c angle*90 
 static void construct_rotation_matrix(int& axis, int& angle,
                                       DenseMatrix<int>& mat);
 
 /// Helper function: Performs the operation : vect2 = mat*vect1
 static void mult_mat_vect(const DenseMatrix<int>& mat,
                           const Vector<int>& vect1,
                           Vector<int>& vect2);
 
 /// Helper function: Performs the operation : mat3=mat1*mat2
 static void mult_mat_mat(const DenseMatrix<int>& mat1,
                          const DenseMatrix<int>& mat2, 
                          DenseMatrix<int>& mat3);


 /// \short Returns the vector of the coordinate directions 
 /// of vertex node number n in an element with nnode1d element per 
 /// dimension.
 static Vector<int> vertex_node_to_vector(const unsigned& n, 
                                          const unsigned& nnode1d);

 

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

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


 /// \short 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)
 static DenseMatrix<bool> Is_adjacent;

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

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

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

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

 /// \short 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]. 
 static DenseMatrix<double> S_steplo; 

 /// \short 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.]
 static DenseMatrix<double> S_stephi; 

 /// \short 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.]
 static DenseMatrix<double> S_directlo;

 /// \short 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.]
 static DenseMatrix<double> S_directhi;

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

 /// \short 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]. 
 static DenseMatrix<double> S_step_edge; 

 /// \short 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). 
 static DenseMatrix<double> S_direct_edge;


};



//===================================================================
/// OcTreeRoot is a OcTree that forms the root of a (recursive)
/// octree. The "root node" is special as it holds additional 
/// information about its neighbours and their relative 
/// rotation (inside a OcTreeForest).
//==================================================================
class OcTreeRoot : public virtual OcTree, public virtual TreeRoot
{
  private:

 /// \short Map of pointers to the edge-neighbouring [Oc]TreeRoots:
 /// Edge_neighbour_pt[direction] is Vector to the pointers to the
 /// [Oc]TreeRoot's edge neighbours in the (enumerated) (edge) direction.
 std::map<int,Vector<TreeRoot*> > Edge_neighbour_pt;

 /// \short Map giving the Up equivalent of the neighbour specified by 
 /// pointer: When viewed from the current octree's neighbour,
 /// our up direction is the neighbour's Up_equivalent[neighbour_pt]
 /// direction. If there's no rotation, this map contains the identify
 /// so that, e.g. \c Up_equivalent[neighbour_pt]=U (read as: "in my
 /// neighbour, my Up is its Up"). If the neighbour is rotated
 /// by 180 degrees relative to the current octree(around the back-front axis),
 /// say, we have \c Up_equivalent[neighbour_pt]=D (read as: "in my 
 /// neighbour, my Up is its Down"); etc.
 std::map<TreeRoot*,int> Up_equivalent;

 /// \short Map giving the Right equivalent of the neighbour specified by 
 /// pointer: When viewed from the current octree's neighbour,
 /// our right direction is the neighbour's right_equivalent[neighbour_pt]
 /// direction. If there's no rotation, this map contains the identify
 /// so that, e.g. \c Right_equivalent[neighbour_pt]=R (read as: "in my
 /// neighbour, my Right is its Right").
 std::map<TreeRoot*,int> Right_equivalent;


  public:

 /// \short Constructor for the root octree: Pass pointer to the 
 /// RefineableQElement<3> that is represented by the OcTree.
 OcTreeRoot(RefineableElement* const &object_pt) : Tree(object_pt),
  OcTree(object_pt), TreeRoot(object_pt)
  {
   
#ifdef PARANOID
   // Check that static member data has been set up
   if (!Static_data_has_been_setup)
    {
    std::string error_message =
     "Static member data hasn't been setup yet.\n";
     error_message +=
      "Call OcTree::setup_static_data() before creating\n";
     error_message += "any OcTreeRoots\n";
 
     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   
  }


 /// Broken copy constructor
 OcTreeRoot(const OcTreeRoot& dummy) : TreeRoot(dummy) 
  { 
   BrokenCopy::broken_copy("OcTreeRoot");
  } 
 
 /// Broken assignment operator
 void operator=(const OcTreeRoot&) 
  {
   BrokenCopy::broken_assign("OcTreeRoot");
  }


 /// \short Return vector of pointers to the edge-neighbouring TreeRoots
 /// in the (enumerated) (edge) direction.
 Vector<TreeRoot*> edge_neighbour_pt(const unsigned& edge_direction)
  {
   
#ifdef PARANOID
   using namespace OcTreeNames;
   if ((edge_direction!=LB)&&(edge_direction!=RB)&&(edge_direction!=DB)&&
       (edge_direction!=UB)&&
       (edge_direction!=LD)&&(edge_direction!=RD)&&(edge_direction!=LU)&&
       (edge_direction!=RU)&&
       (edge_direction!=LF)&&(edge_direction!=RF)&&(edge_direction!=DF)&&
       (edge_direction!=UF))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong edge_direction: " << Direct_string[edge_direction] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   return Edge_neighbour_pt[edge_direction];
  }


 /// \short Return number of edge-neighbouring OcTreeRoot
 /// in the (enumerated) (edge) direction.
 unsigned nedge_neighbour(const unsigned& edge_direction)
 {

#ifdef PARANOID
   using namespace OcTreeNames;
   if ((edge_direction!=LB)&&(edge_direction!=RB)&&(edge_direction!=DB)&&
       (edge_direction!=UB)&&
       (edge_direction!=LD)&&(edge_direction!=RD)&&(edge_direction!=LU)&&
       (edge_direction!=RU)&&
       (edge_direction!=LF)&&(edge_direction!=RF)&&(edge_direction!=DF)&&
       (edge_direction!=UF))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong edge_direction: " << Direct_string[edge_direction] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
  return Edge_neighbour_pt[edge_direction].size();
 }
 
 /// \short Add pointer to the edge-neighbouring [Oc]TreeRoot
 /// in the (enumerated) (edge) direction -- maintains uniqueness
 void add_edge_neighbour_pt(TreeRoot* oc_tree_root_pt,
                            const unsigned& edge_direction)
  {

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

   Vector<TreeRoot*>::iterator it=
    find(Edge_neighbour_pt[edge_direction].begin(),
         Edge_neighbour_pt[edge_direction].end(),
         oc_tree_root_pt);
   if (it==Edge_neighbour_pt[edge_direction].end())
    {
     Edge_neighbour_pt[edge_direction].push_back(oc_tree_root_pt);
    }
  }


 /// \short Return up equivalent of the neighbours specified by
 /// pointer: When viewed from the current octree's neighbour,
 /// our up direction is the neighbour's Up_equivalent[neighbour_pt] 
 /// direction. If there's no rotation, this map contains the identify
 /// so that, e.g. \c Up_equivalent[neighbour_pt]=U (read as: "in my
 /// neighbour, my Up is its Up"). If the neighbour is rotated
 /// by 180 degrees relative to the current octree (around the back-front axis)
 /// say, we have \c Up_equivalent[neighbour_pt]=D (read as: "in my 
 /// neighbour, my Up is its Down"); etc. Returns OMEGA if the Octree 
 /// specified by the pointer argument is not a neighbour. 
 int up_equivalent(TreeRoot* tree_root_pt)
  {
   if (direction_of_neighbour(tree_root_pt)==OMEGA)
    {
     return OMEGA;
    }
   else
    {
     return Up_equivalent[tree_root_pt];
    }
  }


 /// \short Set up equivalent of the neighbours specified by
 /// pointer: When viewed from the current octree's neighbour,
 /// our up direction is the neighbour's Up_equivalent[neighbour_pt] 
 /// direction. If there's no rotation, this map contains the identify
 /// so that, e.g. \c Up_equivalent[neighbour_pt]=U (read as: "in my
 /// neighbour, my Up is its Up"). If the neighbour is rotated
 /// by 180 degrees relative to the current octree (around the back-front axis)
 /// say, we have \c Up_equivalent[neighbour_pt]=D (read as: "in my 
 /// neighbour, my Up is its Down"); etc.
 void set_up_equivalent(TreeRoot* tree_root_pt, const int& dir)
  {
   Up_equivalent[tree_root_pt]=dir;
  }


 /// \short The same thing as up_equivalent, but for the right direction: 
 /// When viewed from the current octree neighbour, our 
 /// right direction is the neighbour's Right_equivalent[neighbour_pt] 
 /// direction. Returns OMEGA if the Octree specified by the pointer
 /// argument is not a neighbour. 
 int right_equivalent(TreeRoot* tree_root_pt)
  {
   if (direction_of_neighbour(tree_root_pt)==OMEGA)
    {
     return OMEGA;
    }
   else
    {
     return Right_equivalent[tree_root_pt];
    }
  }

 /// \short The same thing as up_equivalent, but for the right direction: 
 /// When viewed from the current octree neighbour, our 
 /// right direction is the neighbour's Right_equivalent[neighbour_pt] 
 /// direction.
 void set_right_equivalent(TreeRoot* tree_root_pt, const int& dir)
  {
   Right_equivalent[tree_root_pt]=dir;
  }

 /// \short If octree_root_pt is a neighbour, return the direction
 /// [faces L/R/F/B/U/D or edges DB/UP/...] in which it is found, 
 /// otherwise return OMEGA
 int direction_of_neighbour(TreeRoot* octree_root_pt)
  {
   using namespace OcTreeNames;

   if(Neighbour_pt[U]==octree_root_pt) {return U;}
   if(Neighbour_pt[D]==octree_root_pt) {return D;}
   if(Neighbour_pt[L]==octree_root_pt) {return L;}
   if(Neighbour_pt[R]==octree_root_pt) {return R;}
   if(Neighbour_pt[F]==octree_root_pt) {return F;}
   if(Neighbour_pt[B]==octree_root_pt) {return B;}

   if(Neighbour_pt[LB]==octree_root_pt) {return LB;}
   if(Neighbour_pt[RB]==octree_root_pt) {return RB;}
   if(Neighbour_pt[DB]==octree_root_pt) {return DB;}
   if(Neighbour_pt[UB]==octree_root_pt) {return UB;}

   if(Neighbour_pt[LD]==octree_root_pt) {return LD;}
   if(Neighbour_pt[RD]==octree_root_pt) {return RD;}
   if(Neighbour_pt[LU]==octree_root_pt) {return LU;}
   if(Neighbour_pt[RU]==octree_root_pt) {return RU;}

   if(Neighbour_pt[LF]==octree_root_pt) {return LF;}
   if(Neighbour_pt[RF]==octree_root_pt) {return RF;}
   if(Neighbour_pt[DF]==octree_root_pt) {return DF;}
   if(Neighbour_pt[UF]==octree_root_pt) {return UF;}


   // Search over all edge neighbours
   for (int dir=LB;dir<=UF;dir++)
    {
     Vector<TreeRoot*> edge_neigh_pt=this->edge_neighbour_pt(dir);
     unsigned n_neigh=edge_neigh_pt.size();
     for (unsigned e=0;e<n_neigh;e++)
      {
       if (edge_neigh_pt[e]==octree_root_pt) 
        {
         return dir;
        }
      }
    }


   //If we get here, it's not a neighbour
   return OMEGA;
  }

};
 


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

 
//================================================================
/// An OcTreeForest consists of a collection of OcTreeRoots.
/// Each member tree can have neighbours to its L/R/U/D/F/B and
/// DB/UP/... and the orientation of their compasses can differ, 
/// allowing for complex, unstructured meshes. 
//=================================================================
class OcTreeForest : public TreeForest
{ 

 public:

 /// \short Constructor for OcTree forest: Pass Vector of 
 /// (pointers to) trees.
 OcTreeForest(Vector<TreeRoot* >& trees_pt);

 /// Default constructor (empty and broken)
 OcTreeForest()
  {
   throw OomphLibError("Don't call empty constructor for OcTreeForest!",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

 /// Broken copy constructor
 OcTreeForest(const OcTreeForest& dummy) 
  { 
   BrokenCopy::broken_copy("OcTreeForest");
  } 
 
 /// Broken assignment operator
 void operator=(const OcTreeForest&) 
  {
   BrokenCopy::broken_assign("OcTreeForest");
  }

 /// \short Destructor: Delete the constituent octrees (and thus 
 /// the associated objects!)
 virtual ~OcTreeForest() {}
 


 /// \short 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.disable_doc() 
 /// has been called, no output is created.
 void check_all_neighbours(DocInfo &doc_info);
 
/// \short Open output files that will store any hanging nodes in
 /// the forest and return a vector of the streams.
 void open_hanging_node_files(DocInfo &doc_info,
                              Vector<std::ofstream*> &output_stream);

 /// \short Self-test: Check all neighbours. Return success (0) 
 /// if the max. distance between corresponding points in the
 /// neighbours is less than the tolerance specified in the
 /// static value Tree::Max_neighbour_finding_tolerance.
 unsigned self_test();

 /// \short Return pointer to i-th OcTree in forest
 /// (Performs a dynamic cast from the TreeRoot to a 
 /// OcTreeRoot). 
 OcTreeRoot* octree_pt(const unsigned &i) const 
  {return dynamic_cast<OcTreeRoot*>(Trees_pt[i]);}
 
 /// \short Given the number i of the root octree in this forest, return
 /// pointer to its face neighbour in the specified direction. NULL
 /// if neighbour doesn't exist. (This does the dynamic cast 
 /// from a TreeRoot to a OcTreeRoot internally).
 OcTreeRoot* oc_face_neigh_pt(const unsigned &i, const int &direction)
  {
#ifdef PARANOID
   using namespace OcTreeNames;
   if ((direction!=U)&&(direction!=D)&&
       (direction!=F)&&(direction!=B)&&
       (direction!=L)&&(direction!=R))
    {
     std::ostringstream error_stream;
     error_stream << "Wrong edge_direction: " 
                  << OcTree::Direct_string[direction] 
                  << std::endl;
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return dynamic_cast<OcTreeRoot*>(Trees_pt[i]->neighbour_pt(direction));
  }

 /// \short Given the number i of the root octree in this forest, return
 /// the vector of pointers to the true edge neighbours in the specified 
 /// (edge) direction.
 Vector<TreeRoot*> oc_edge_neigh_pt(const unsigned &i, const int &direction)
  {
   // Note: paranoia check is done in edge_neighbour_pt
   return dynamic_cast<OcTreeRoot*>(
    Trees_pt[i])->edge_neighbour_pt(direction);
  }
 

private:


 /// Construct the rotation schemes
 void construct_up_right_equivalents();
 
 /// Construct the neighbour scheme
 void find_neighbours();
 

};

}

#endif