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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#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
  {
    /// 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
    };
  } // namespace OcTreeNames
 
 
  // 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:
    /// Destructor. Note: Deleting a octree also deletes the
    /// objects associated with all non-leaf nodes!
    virtual ~OcTree() {}
 
 
    /// Broken copy constructor
    OcTree(const OcTree& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const OcTree&) = delete;
 
    /// 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;
    }
 
    /// Function that, given an edge, returns the two faces on which it
    // lies between, i.e. the faces to which it is a common edge
    static Vector<int> faces_of_common_edge(const int& edge);
 
    /// 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,
                                bool& in_neighbouring_tree) const;
 
    /// 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;
 
 
    /// 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();
 
    /// Setup the static data, rotation and reflection schemes, etc
    static void setup_static_data();
 
 
    /// 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);
 
 
    /// 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);
 
    /// 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);
 
    /// 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);
 
 
    /// 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);
 
 
    /// 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);
 
 
    /// 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;
 
    /// \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;
 
    /// 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;
 
    /// 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;
 
 
    /// 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);
    }
 
    /// 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) {}
 
 
    /// 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:
    /// 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,
                                bool& in_neighbouring_tree,
                                int max_level,
                                OcTreeRoot* orig_root_pt) const;
 
 
    /// 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;
 
 
    /// 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;
 
 
    /// 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);
 
 
    /// 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;
 
 
    /// 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;
 
    /// Reflection scheme: Reflect(direction,octant): Get mirror
    /// of octant/edge in specified direction. E.g. Reflect(LDF,L)=RDF
    static DenseMatrix<int> Reflect;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// S_base_edge(i,edge):  Initial value for coordinate s[i] on
    /// the specified edge (LF/RF/...).
    static DenseMatrix<double> 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].
    static DenseMatrix<double> 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).
    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:
    /// 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;
 
    /// 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;
 
    /// 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:
    /// 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) = delete;
 
    /// Broken assignment operator
    void operator=(const OcTreeRoot&) = delete;
 
    /// 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];
    }
 
 
    /// 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();
    }
 
    /// 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);
      }
    }
 
 
    /// 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];
      }
    }
 
 
    /// 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;
    }
 
 
    /// 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];
      }
    }
 
    /// 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;
    }
 
    /// 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:
    /// 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) = delete;
 
    /// Broken assignment operator
    void operator=(const OcTreeForest&) = delete;
 
    /// Destructor: Delete the constituent octrees (and thus
    /// the associated objects!)
    virtual ~OcTreeForest() {}
 
 
    /// 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);
 
    /// 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);
 
    /// 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();
 
    /// 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]);
    }
 
    /// 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));
    }
 
    /// 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);
    }
 
    /// Construct the rotation schemes
    void construct_up_right_equivalents();
 
  private:
    /// Construct the neighbour scheme
    void find_neighbours();
  };
 
} // namespace oomph
 
#endif