binary_tree.cc

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Non-inline and non-templated functions for BinaryTree and BinaryTreeForest
// classes
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include <set>
 
// oomph-lib headers
#include "binary_tree.h"
#include "refineable_line_element.h"
 
namespace oomph
{
  //========================================================================
  /// Boolean indicating that static member data has been setup.
  //========================================================================
  bool BinaryTree::Static_data_has_been_setup = false;
 
  //========================================================================
  /// Colours for neighbours in various directions (static data).
  //========================================================================
  Vector<std::string> BinaryTree::Colour;
 
  //========================================================================
  /// S_base(direction): Initial value for coordinate s on the edge
  /// indicated by direction (L/R) (static data).
  //========================================================================
  Vector<double> BinaryTree::S_base;
 
  //========================================================================
  /// Translate (enumerated) directions into strings (static data).
  //========================================================================
  Vector<std::string> BinaryTree::Direct_string;
 
  //========================================================================
  /// Get opposite edge, e.g. Reflect_edge[N]=S (static data)
  //========================================================================
  Vector<int> BinaryTree::Reflect_edge;
 
  //========================================================================
  /// Array of direction/segment adjacency scheme:
  /// Is_adjacent(i_vertex,j_segment): Is vertex adjacent to segment?
  /// (static data)
  //========================================================================
  DenseMatrix<bool> BinaryTree::Is_adjacent;
 
  //========================================================================
  /// Reflection scheme: Reflect(direction,segment): Get mirror of segment
  /// in specified direction. E.g. Reflect(L,L)=R (static data)
  //========================================================================
  DenseMatrix<int> BinaryTree::Reflect;
 
  //========================================================================
  /// Set up the static data stored in the BinaryTree -- this needs to be
  /// called before BinaryTrees can be used. Automatically called by
  /// RefineableLineMesh constructor.
  //========================================================================
  void BinaryTree::setup_static_data()
  {
    using namespace BinaryTreeNames;
 
#ifdef PARANOID
    if (Tree::OMEGA != BinaryTree::OMEGA)
    {
      std::ostringstream error_stream;
      error_stream << "Inconsistent enumeration!  \n    Tree::OMEGA="
                   << Tree::OMEGA << "\nBinaryTree::OMEGA=" << BinaryTree::OMEGA
                   << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Set flag to indicate that static data has been setup
    Static_data_has_been_setup = true;
 
    // Tecplot colours for neighbours in various directions
    Colour.resize(27);
 
    Colour[L] = "RED";
    Colour[R] = "CYAN";
    Colour[OMEGA] = "YELLOW";
 
    // S_base(direction): Initial value for coordinate s on the
    // edge indicated by direction (L/R)
    S_base.resize(27);
 
    S_base[L] = -1.0;
    S_base[R] = 1.0;
 
    // Translate (enumerated) directions into strings
    Direct_string.resize(27);
 
    Direct_string[L] = "L";
    Direct_string[R] = "R";
    Direct_string[OMEGA] = "OMEGA";
 
    // Build direction/segment adjacency scheme:
    // Is_adjacent(i_vertex,j_segment): Is vertex adjacent to segment?
    Is_adjacent.resize(27, 27);
 
    Is_adjacent(L, L) = true;
    Is_adjacent(R, L) = false;
    Is_adjacent(L, R) = false;
    Is_adjacent(R, R) = true;
 
    // Build reflection scheme: Reflect(direction,segment):
    // Get mirror of segment in direction
    Reflect.resize(27, 27);
 
    Reflect(L, L) = R;
    Reflect(R, L) = R;
    Reflect(L, R) = L;
    Reflect(R, R) = L;
 
    // Get opposite edge, e.g. Reflect_edge(L)=R
    Reflect_edge.resize(27);
 
    Reflect_edge[L] = R;
    Reflect_edge[R] = L;
  }
 
  //========================================================================
  /// Return pointer to greater or equal-sized edge neighbour in specified
  /// \c direction; also provide info regarding the relative size of the
  /// neighbour:
  /// - In the present binary tree, the left vertex is located at the local
  ///   coordinate s = -1. This point is located at the local coordinate
  ///   s = \c s_in_neighbour[0] in the neighbouring binary tree.
  /// - We're looking for a neighbour in the specified \c direction. When
  ///   viewed from the neighbouring binary tree, the edge that separates
  ///   the present binary tree from its neighbour is the neighbour's
  ///   \c edge edge. Since in 1D there can be no rotation between the two
  ///   binary trees, this is a simple reflection. For instance, if we're
  ///   looking for a neighhbour in the \c L [eft] \c direction, \c edge
  ///   will be \c R [ight].
  /// - \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 binary tree.
  /// - \c in_neighbouring_tree returns true if we have had to flip to a
  ///   different root, even if that root is actually the same (as it can
  ///   be in periodic problems).
  //========================================================================
  BinaryTree* BinaryTree::gteq_edge_neighbour(const int& direction,
                                              Vector<double>& s_in_neighbour,
                                              int& edge,
                                              int& diff_level,
                                              bool& in_neighbouring_tree) const
  {
    using namespace BinaryTreeNames;
 
#ifdef PARANOID
    if ((direction != L) && (direction != R))
    {
      std::ostringstream error_stream;
      error_stream << "Direction " << direction << " is not L or R"
                   << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Initialise in_neighbouring tree to false. It will be set to true
    // during the recursion if we do actually hop over into the neighbour
    in_neighbouring_tree = false;
 
    // Maximum level to which we're allowed to descend (we only want
    // greater-or-equal-sized neighbours)
    int max_level = Level;
 
    // Current element has the following root:
    BinaryTreeRoot* orig_root_pt = dynamic_cast<BinaryTreeRoot*>(Root_pt);
 
    // Initialise offset in local coordinate
    double s_diff = 0;
 
    // Initialise difference in level
    diff_level = 0;
 
    // Find neighbour
    BinaryTree* return_pt = gteq_edge_neighbour(direction,
                                                s_diff,
                                                diff_level,
                                                in_neighbouring_tree,
                                                max_level,
                                                orig_root_pt);
 
    BinaryTree* neighb_pt = return_pt;
 
    // If neighbour exists...
    if (neighb_pt != 0)
    {
      // What's the direction of the interfacial edge when viewed from within
      // the neighbour element?
      edge = Reflect_edge[direction];
 
      // What's the local coordinate s at which this edge is located in the
      // neighbouring element?
      s_in_neighbour[0] = S_base[edge];
    }
    return return_pt;
  }
 
  //========================================================================
  /// Find `greater-or-equal-sized edge neighbour' in given direction (L/R).
  ///
  /// This is an auxiliary routine which allows neighbour finding in
  /// adjacent binary trees. Needs to keep track of previous son types
  /// and the maximum level to which search is performed.
  ///
  /// Parameters:
  /// - direction (L/R): Direction in which neighbour has to be found.
  /// - s_diff: Offset of left vertex from corresponding vertex in
  ///   neighbour. Note that this is input/output as it needs to be
  ///   incremented/decremented during the recursive calls to this function.
  /// - edge: We're looking for the neighbour across our edge 'direction'
  ///   (L/R). When viewed from the neighbour, this edge is `edge' (L/R).
  ///   Since there is no relative rotation between neighbours this is a
  ///   mere reflection, e.g. direction=L  --> edge=R etc.
  /// - diff_level <= 0 indicates the difference in binary tree levels
  ///   between the current element and its neighbour.
  /// - max_level is the maximum level to which the neighbour search is
  ///   allowed to proceed. This is again necessary because in a forest,
  ///   the neighbour search isn't based on pure recursion.
  /// - orig_root_pt identifies the root node of the element whose
  ///   neighbour we're really trying to find by all these recursive calls.
  //========================================================================
  BinaryTree* BinaryTree::gteq_edge_neighbour(
    const int& direction,
    double& s_diff,
    int& diff_level,
    bool& in_neighbouring_tree,
    int max_level,
    BinaryTreeRoot* const& orig_root_pt) const
  {
    using namespace BinaryTreeNames;
 
#ifdef PARANOID
    if ((direction != L) && (direction != R))
    {
      std::ostringstream error_stream;
      error_stream << "Direction " << direction << " is not L or R"
                   << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    BinaryTree* next_el_pt;
    BinaryTree* return_el_pt;
 
    // STEP 1: Find the neighbour's father
    // -------
 
    // Does the element have a father?
    if (Father_pt != 0)
    {
      // If the present segment (whose location inside its father element
      // is specified by Son_type) is adjacent to the father's edge in the
      // required direction, then its neighbour has a different father
      // ---> we need to climb up the tree to the father and find his
      // neighbour in the required direction. Note that this is the cunning
      // recursive part. The returning may not stop until we hit the very
      // top of the tree, when the element does NOT have a father.
      if (Is_adjacent(direction, Son_type))
      {
        next_el_pt = dynamic_cast<BinaryTree*>(Father_pt)->gteq_edge_neighbour(
          direction,
          s_diff,
          diff_level,
          in_neighbouring_tree,
          max_level,
          orig_root_pt);
      }
      // If the present segment is not adjacent to the father's edge in the
      // required direction, then the neighbour has the same father and is
      // obtained by the appropriate reflection inside the father element.
      // This will only be called if we have not left the original tree.
      else
      {
        next_el_pt = dynamic_cast<BinaryTree*>(Father_pt);
      }
 
      // We're about to ascend one level
      diff_level -= 1;
 
      // Work out position of left corner of present edge in its father element
      s_diff += pow(0.5, -diff_level);
 
      // STEP 2: We have now located the neighbour's father and need to
      // ------- find the appropriate son.
 
      // Buffer cases where the neighbour (and hence its father) lie outside
      // the boundary
      if (next_el_pt != 0)
      {
        // If the father is a leaf then we can't descend to the same
        // level as the present node ---> simply return the father himself
        // as the (greater) neighbour. Same applies if we are about to
        // descend lower than the max_level (in a  neighbouring tree).
        if ((next_el_pt->Son_pt.size() == 0) ||
            (next_el_pt->Level > max_level - 1))
        {
          return_el_pt = next_el_pt;
        }
        // We have located the neighbour's father: The position of the
        // neighbour is obtained by `reflecting' the position of the
        // node itself. We know exactly how to reflect because we know which
        // son type we are and we have the pointer to the neighbour's father.
        else
        {
          int son_segment = Reflect(direction, Son_type);
 
          // The next element in the tree is the appropriate son of the
          // neighbour's father
          return_el_pt =
            dynamic_cast<BinaryTree*>(next_el_pt->Son_pt[son_segment]);
 
          // Work out position of left corner of present edge in next
          // higher element
          s_diff -= pow(0.5, -diff_level);
 
          // We have just descended one level
          diff_level += 1;
        }
      }
      // The neighbour's father lies outside the boundary --> the neighbour
      // itself does too --> return NULL
      else
      {
        return_el_pt = 0;
      }
    }
    // Element does not have a father --> check if it has a neighbouring
    // tree in the appropriate direction
    else
    {
      // Find neighbouring root
      if (Root_pt->neighbour_pt(direction) != 0)
      {
        // In this case we have moved to a neighbour, so set the flag
        in_neighbouring_tree = true;
        return_el_pt =
          dynamic_cast<BinaryTreeRoot*>(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;
  }
 
  //========================================================================
  /// 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 neighbour. If the difference is less than
  /// Tree::Max_neighbour_finding_tolerance return success (0), otherwise
  /// failure (1).
  //========================================================================
  unsigned BinaryTree::self_test()
  {
    // Stick pointers to all nodes into Vector and number elements
    // in the process
    Vector<Tree*> all_nodes_pt;
    stick_all_tree_nodes_into_vector(all_nodes_pt);
    long int count = 0;
    unsigned long num_nodes = all_nodes_pt.size();
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      all_nodes_pt[i]->object_pt()->set_number(++count);
    }
 
    // Check neighbours (distance between hanging nodes) -- don't print
    // (keep output streams closed)
    double max_error = 0.0;
    std::ofstream neighbours_file;
    std::ofstream neighbours_txt_file;
    BinaryTree::doc_neighbours(
      all_nodes_pt, neighbours_file, neighbours_txt_file, max_error);
 
    if (max_error > Max_neighbour_finding_tolerance)
    {
      oomph_info << "\n \n Failed self_test() for BinaryTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 1;
    }
    else
    {
      oomph_info << "\n \n Passed self_test() for BinaryTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 0;
    }
  }
 
 
  //========================================================================
  /// Constructor for BinaryTreeForest:
  ///
  /// Pass:
  ///  - trees_pt[], the Vector of pointers to the constituent trees
  ///    (BinaryTreeRoot objects)
  //========================================================================
  BinaryTreeForest::BinaryTreeForest(Vector<TreeRoot*>& trees_pt)
    : TreeForest(trees_pt)
  {
#ifdef LEAK_CHECK
    LeakCheckNames::BinaryTreeForest_build += 1;
#endif
 
    using namespace BinaryTreeNames;
 
    // Set up the neighbours
    find_neighbours();
  }
 
  //========================================================================
  /// Set up the neighbour lookup schemes for all constituent binary trees.
  //========================================================================
  void BinaryTreeForest::find_neighbours()
  {
    using namespace BinaryTreeNames;
 
    unsigned numtrees = ntree();
    unsigned n = 0; // to store nnode1d
    if (numtrees > 0)
    {
      n = Trees_pt[0]->object_pt()->nnode_1d();
    }
    else
    {
      throw OomphLibError(
        "Trying to setup the neighbour scheme for an empty forest\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    // Number of vertex nodes: 2
    unsigned n_vertex_node = 2;
 
    // Find 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<LineElementBase*>(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 neighbouring trees
    // i.e. trees that share a node
    Vector<std::set<unsigned>> 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 connected
          if (i != j)
          {
            neighb_tree[i].insert(j);
          }
        }
      }
    }
 
    // Loop over all trees
    for (unsigned i = 0; i < numtrees; i++)
    {
      // Loop over their neighbours
      for (std::set<unsigned>::iterator it = neighb_tree[i].begin();
           it != neighb_tree[i].end();
           it++)
      {
        unsigned j = (*it);
 
        // is it the left-hand neighbour?
        bool is_L_neighbour = Trees_pt[j]->object_pt()->get_node_number(
                                Trees_pt[i]->object_pt()->node_pt(0)) != -1;
 
        // is it the right-hand neighbour?
        bool is_R_neighbour = Trees_pt[j]->object_pt()->get_node_number(
                                Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1;
 
        if (is_L_neighbour) Trees_pt[i]->neighbour_pt(L) = Trees_pt[j];
        if (is_R_neighbour) Trees_pt[i]->neighbour_pt(R) = Trees_pt[j];
      }
    } // End of loop over all trees
  }
 
  //========================================================================
  /// Document and check all the neighbours in all the nodes in the forest.
  //========================================================================
  void BinaryTreeForest::check_all_neighbours(DocInfo& doc_info)
  {
    // Create Vector of elements
    Vector<Tree*> all_tree_nodes_pt;
    this->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
 
    // Create storage for information files
    std::ofstream neigh_file;
    std::ofstream neigh_txt_file;
 
    // If we are documenting the results, then open the files
    if (doc_info.is_doc_enabled())
    {
      std::ostringstream fullname;
      fullname << doc_info.directory() << "/neighbours" << doc_info.number()
               << ".dat";
      oomph_info << "opened " << fullname.str() << " to doc neighbours"
                 << std::endl;
      neigh_file.open(fullname.str().c_str());
      fullname.str("");
      fullname << doc_info.directory() << "/neighbours" << doc_info.number()
               << ".txt";
      oomph_info << "opened " << fullname.str() << " to doc neighbours"
                 << std::endl;
      neigh_txt_file.open(fullname.str().c_str());
    }
 
    // Call the standard documentation function
    double max_error = 0.0;
    BinaryTree::doc_neighbours(
      all_tree_nodes_pt, neigh_file, neigh_txt_file, max_error);
 
    // If the error is too large, complain
    if (max_error > Tree::max_neighbour_finding_tolerance())
    {
      std::ostringstream error_stream;
      error_stream << "Max. error in binary tree neighbour finding: "
                   << max_error << " is too big" << std::endl;
      error_stream
        << "i.e. bigger than Tree::max_neighbour_finding_tolerance()="
        << Tree::max_neighbour_finding_tolerance() << std::endl;
 
      // Close the files if they were opened
      if (doc_info.is_doc_enabled())
      {
        neigh_file.close();
        neigh_txt_file.close();
      }
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    else
    {
      oomph_info << "Max. error in binary tree neighbour finding: " << max_error
                 << " is OK" << std::endl;
      oomph_info
        << "i.e. less than BinaryTree::max_neighbour_finding_tolerance()="
        << BinaryTree::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();
    }
  }
 
  //========================================================================
  /// 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 neighbour. If the difference is less than
  /// Tree::Max_neighbour_finding_tolerance return success (0), otherwise
  /// failure (1).
  //========================================================================
  unsigned BinaryTreeForest::self_test()
  {
    // Stick pointers to all nodes into Vector and number elements
    // in the process
    Vector<Tree*> all_forest_nodes_pt;
    stick_all_tree_nodes_into_vector(all_forest_nodes_pt);
    long int count = 0;
    unsigned long num_nodes = all_forest_nodes_pt.size();
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      all_forest_nodes_pt[i]->object_pt()->set_number(++count);
    }
 
    // Check neighbours (distance between hanging nodes) -- don't print
    // (keep output streams closed)
    double max_error = 0.0;
    std::ofstream neighbours_file;
    std::ofstream neighbours_txt_file;
    BinaryTree::doc_neighbours(
      all_forest_nodes_pt, neighbours_file, neighbours_txt_file, max_error);
    if (max_error > BinaryTree::max_neighbour_finding_tolerance())
    {
      oomph_info << "\n \n Failed self_test() for BinaryTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 1;
    }
    else
    {
      oomph_info << "\n \n Passed self_test() for BinaryTree: Max. error "
                 << max_error << std::endl
                 << std::endl;
      return 0;
    }
  }
 
  //========================================================================
  /// Doc/check all neighbours of binary tree ("nodes") contained in the
  /// Vector forest_node_pt. Output into neighbours_file which can be viewed
  /// from tecplot with BinaryTreeNeighbours.mcr. Neighbour info and errors
  /// are displayed on neighbours_txt_file. Finally, compute the maximum
  /// error between vertices when viewed from neighbouring element. Output
  /// is suppressed if the output streams are closed.
  //========================================================================
  void BinaryTree::doc_neighbours(Vector<Tree*> forest_nodes_pt,
                                  std::ofstream& neighbours_file,
                                  std::ofstream& neighbours_txt_file,
                                  double& max_error)
  {
    using namespace BinaryTreeNames;
 
    int diff_level;
    double s_diff;
    bool in_neighbouring_tree;
    int edge = OMEGA;
 
    Vector<double> s(1);
    Vector<double> x(1);
    Vector<int> prev_son_type;
 
    Vector<double> s_in_neighbour(1);
 
    Vector<double> x_small(1);
    Vector<double> x_large(1);
 
    // Initialise error in vertex positions
    max_error = 0.0;
 
    // Loop over all elements to assign numbers for plotting
    // -----------------------------------------------------
    unsigned long num_nodes = forest_nodes_pt.size();
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      // Set number
      forest_nodes_pt[i]->object_pt()->set_number(i);
    }
 
    // Loop over all elements for checks
    // ---------------------------------
    for (unsigned long i = 0; i < num_nodes; i++)
    {
      // Doc the element itself
      BinaryTree* el_pt = dynamic_cast<BinaryTree*>(forest_nodes_pt[i]);
 
      // If the object is incomplete, complain
      if (el_pt->object_pt()->nodes_built())
      {
        // Print it
        if (neighbours_file.is_open())
        {
          dynamic_cast<RefineableQElement<1>*>(el_pt->object_pt())
            ->output_corners(neighbours_file, "BLACK");
        }
 
        // Loop over directions to find neighbours
        // ---------------------------------------
        for (int direction = L; direction <= R; direction++)
        {
          // Initialise difference in levels and coordinate offset
          diff_level = 0;
          s_diff = 0.0;
 
          // Find greater-or-equal-sized neighbour...
          BinaryTree* neighb_pt = el_pt->gteq_edge_neighbour(
            direction, s_in_neighbour, edge, diff_level, in_neighbouring_tree);
 
          // If neighbour exist and nodes are created: Doc it
          if ((neighb_pt != 0) && (neighb_pt->object_pt()->nodes_built()))
          {
            // Doc neighbour stats
            if (neighbours_txt_file.is_open())
            {
              neighbours_txt_file
                << Direct_string[direction] << " neighbour of "
                << el_pt->object_pt()->number() << " is "
                << neighb_pt->object_pt()->number() << " diff_level "
                << diff_level << " s_diff " << s_diff
                << " inside neighbour the edge is " << Direct_string[edge]
                << std::endl
                << std::endl;
            }
 
            // Plot neighbour in the appropriate colour
            if (neighbours_file.is_open())
            {
              dynamic_cast<RefineableQElement<1>*>(neighb_pt->object_pt())
                ->output_corners(neighbours_file, Colour[direction]);
            }
 
            // Check that local coordinates in the larger element (obtained
            // via s_diff) lead to the same spatial point as the node vertices
            // in the current element
            {
              if (neighbours_file.is_open())
              {
                neighbours_file << "ZONE I=1 \n";
              }
 
              // Left vertex:
              // ------------
 
              // Get coordinates in large (neighbour) element
              s[0] = s_in_neighbour[0];
              neighb_pt->object_pt()->get_x(s, x_large);
 
              // Get coordinates in small element
              Vector<double> s(1);
              s[0] = S_base[direction];
              el_pt->object_pt()->get_x(s, x_small);
 
              // Need to exclude periodic nodes from this check. There can
              // only be periodic nodes if we have moved into the neighbour
              bool is_periodic = false;
              if (in_neighbouring_tree)
              {
                // Is the node periodic?
                is_periodic =
                  el_pt->root_pt()->is_neighbour_periodic(direction);
              }
 
              double error = 0.0;
              // Only bother to calculate the error if the node is NOT periodic
              if (is_periodic == false)
              {
                error += pow(x_small[0] - x_large[0], 2);
              }
 
              // Take the root of the square error
              error = sqrt(error);
              if (neighbours_txt_file.is_open())
              {
                neighbours_txt_file << "Error (1) " << error << std::endl;
              }
 
              if (std::fabs(error) > max_error)
              {
                max_error = std::fabs(error);
              }
 
              if (neighbours_file.is_open())
              {
                neighbours_file << x_large[0] << "  0 \n";
              }
 
              // Right vertex:
              // -------------
 
              // Get coordinates in large (neighbour) element
              s[0] = s_in_neighbour[0];
              neighb_pt->object_pt()->get_x(s, x_large);
 
              // Get coordinates in small element
              s[0] = S_base[direction];
              el_pt->object_pt()->get_x(s, x_small);
 
              error = 0.0;
              // Only do this if we are NOT periodic
              if (is_periodic == false)
              {
                error += pow(x_small[0] - x_large[0], 2);
              }
              // Take the root of the square error
              error = sqrt(error);
 
              // error =
              // sqrt(pow(x_small[0]-x_large[0],2)+pow(x_small[1]-x_large[1],2));
              if (neighbours_txt_file.is_open())
              {
                neighbours_txt_file << "Error (2) " << error << std::endl;
              }
 
              if (std::fabs(error) > max_error)
              {
                max_error = std::fabs(error);
              }
 
              if (neighbours_file.is_open())
              {
                neighbours_file << x_large[0] << "  0 \n";
              }
            }
            //        else
            //         {
            //          // No neighbour: write dummy zone so tecplot can find
            //          four
            //          // neighbours for every element
            //          if (neighbours_file.is_open())
            //            {
            //             neighbours_file << "ZONE I=1 \n";
            //             neighbours_file << "-0.05 -0.05   0 \n";
            //             neighbours_file << "-0.05 -0.05   0 \n";
            //            }
            //         }
          }
          // If neighbour does not exist: Insert blank zones into file
          // so that tecplot can find four neighbours for every element
          else
          {
            if (neighbours_file.is_open())
            {
              neighbours_file << "ZONE \n 0.00  0 \n";
              neighbours_file << "ZONE I=1 \n";
              neighbours_file << "-0.05  0 \n";
              neighbours_file << "-0.05  0 \n";
            }
          }
        }
      } // End of case when element can be documented
    }
  }
 
} // namespace oomph