quadtree.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//====================================================================
//Header file for quadtree and quadtree forest classes
#ifndef OOMPH_QUADTREE_HEADER
#define OOMPH_QUADTREE_HEADER


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

#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

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



namespace oomph
{

//====================================================================
/// Namespace for QuadTree directions
//====================================================================
namespace QuadTreeNames
{
 /// \short Directions. OMEGA is used if a direction is undefined 
 /// in a certain context
 enum{SW,SE,NW,NE,N,E,S,W,OMEGA=26};
};

// Forward class definition for class representing the root of a QuadTree
class QuadTreeRoot;

//================================================================
/// QuadTree class: Recursively defined, generalised quadtree.
///
/// A QuadTree has:
/// - a pointer to the object (of type RefineableQElement<2>) that it 
///   represents in a mesh refinement context.
/// - Vector of pointers to its four (SW/SE/NW/NE) sons (which are 
///   themselves quadtrees). 
///   If the Vector of pointers to the sons has zero length,
///   the QuadTree is a "leaf node" in the overall quadtree.
/// - a pointer to its father. If this pointer is NULL, the QuadTree is the
///   the root node of the overall quadtree.
/// 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 quadtrees.
/// 
/// The objects contained in the quadtree are assumed to be
/// (topologically) rectangular elements whose geometry is
/// parametrised by local coordinates \f$ {\bf s} \in [-1,1]^2 \f$.
///
/// The tree can be traversed and actions performed
/// at all its "nodes" or only at the leaf "nodes" ("nodes" without sons).
///
/// Finally, the leaf "nodes" can be split depending on
/// a criteria defined by the object.
///
/// Note that QuadTrees are only generated by splitting existing
/// QuadTrees. Therefore, the constructors are protected. The
/// only QuadTree that "Joe User" can create is
/// the (derived) class QuadTreeRoot.
//=================================================================
class QuadTree : public virtual Tree
{ 
  public:
 
 /// \short Destructor. Note: Deleting a quadtree also deletes the 
 /// objects associated with all non-leaf nodes!
 virtual ~QuadTree() {}

 /// Broken copy constructor
 QuadTree(const QuadTree& dummy) 
  { 
   BrokenCopy::broken_copy("QuadTree");
  } 
 
 /// Broken assignment operator
 void operator=(const QuadTree&) 
  {
   BrokenCopy::broken_assign("QuadTree");
  }
 
 /// \short Overload the function construct_son to ensure that the son
 /// is a specific QuadTree and not a general Tree.
 Tree* construct_son(RefineableElement* const &object_pt,
                     Tree* const &father_pt, const int &son_type)
  {
   QuadTree* temp_quad_pt = new QuadTree(object_pt,father_pt,son_type);
   return temp_quad_pt;
  }

 /// \short Return pointer to greater or equal-sized edge neighbour 
 /// in specified \c direction; also provide info regarding the relative 
 /// size and orientation of neighbour:
 /// - The vector translate_s turns the index of the local coordinate
 ///   in the present quadtree into that of the neighbour. If there are no
 ///   rotations then translate_s[i] = i, but if, for example, the neighbour's
 ///   eastern face is adjacent to our northern face translate_s[0] = 1
 ///   and translate_s[1] = 0. Of course, this could be deduced after the
 ///   fact, but it's easier to do it here.
 /// - In the present quadtree, the lower left (south west) vertex is
 ///   located at local coordinates (-1,-1). This point is located
 ///   at the local coordinates (\c s_lo[0], \c s_lo[1]) in the neighbouring
 ///   quadtree.
 /// - ditto with s_hi: In the present quadtree, the upper right (north east) 
 ///   vertex is located at local coordinates (1,1). This point is located
 ///   at the local coordinates (\c s_hi[0], \c s_hi[1]) in the neighbouring
 ///   quadtree.
 /// - We're looking for a neighbour in the specified \c direction. When
 ///   viewed from the neighbouring quadtree, the edge that separates 
 ///   the present quadtree from its neighbour is the neighbour's \c edge 
 ///   edge. If there's no rotation between the two quadtrees, this is a 
 ///   simple reflection: For instance, if we're looking
 ///   for a neighhbour in the \c N [orthern] \c direction, \c edge will 
 ///   be \c S [outh]
 /// - \c diff_level <= 0 indicates the difference in refinement levels between
 ///   the two neighbours. If \c diff_level==0, the neighbour has the
 ///   same size as the current quadtree.
 /// - \c in_neighbouring_tree indicates whether the neighbour is actually
 ///   in another tree in the forest. The introduction of this flag 
 ///   was necessitated by periodic problems where a TreeRoot can be its
 ///   own neighbour.
 QuadTree* gteq_edge_neighbour(const int& direction, 
                               Vector<unsigned> &translate_s, 
                               Vector<double>& s_lo,  Vector<double>& s_hi, 
                               int& edge, int& diff_level,
                               bool &in_neighbouring_tree) const;
 
 /// \short Traverse Tree: Preorder traverse and stick pointers to
 /// neighbouring leaf nodes (only) into Vector
 void stick_neighbouring_leaves_into_vector(
       Vector<const QuadTree*>& tree_neighbouring_nodes,
       Vector<Vector<double> >& tree_neighbouring_s_lo,
       Vector<Vector<double> >& tree_neighbouring_s_hi,
       Vector<int>& tree_neighbouring_diff_level,
       const QuadTree* my_neigh_pt,
       const int& direction) const;

 /// \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 QuadTree::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 neighbours of quadtree (nodes) contained in the
 /// Vector forest_node_pt. Output into neighbours_file which can
 /// be viewed from tecplot with QuadTreeNeighbours.mcr
 /// Neighbour info and errors are displayed on 
 /// neighbours_txt_file.  Finally, compute the max. error between 
 /// vertices when viewed from neighhbouring element. 
 /// If the two filestreams are closed, output is suppressed.
 static void doc_neighbours(Vector<Tree*> forest_nodes_pt, 
                            std::ofstream& neighbours_file, 
                            std::ofstream& neighbours_txt_file,
                            double& max_error);


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

  protected:
 

 /// Default constructor (empty and broken)
 QuadTree()
  {
   throw OomphLibError("Don't call an empty constructor for a QuadTree object",
                       OOMPH_CURRENT_FUNCTION,OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Default constructor for empty (root) tree: 
 /// no father, no sons; just pass a pointer to its object
 /// Protected because QuadTrees can only be created internally,
 /// during the split operation. Only QuadTreeRoots can be
 /// created externally. 
 QuadTree(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 (SE/SW/NE/NW) it is.
 /// Protected because QuadTrees can only be created internally,
 /// during the split operation.  Only QuadTreeRoots can be
 /// created externally. 
 QuadTree(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 edge neighbour in direction.
 /// Auxiliary internal routine which passes additional information around.
 QuadTree* gteq_edge_neighbour(const int& direction, 
                               double& s_diff, 
                               int& diff_level,
                               bool& in_neighbouring_tree,
                               int max_level, 
                               QuadTreeRoot* const &orig_root_pt) const;

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

 /// \short S_base(i,direction):  Initial value for coordinate s[i] on 
 /// the edge indicated by direction (S/E/N/W) 
 static DenseMatrix<double> S_base;
 
 /// \short S_step(i,direction) Increments for coordinate s[i] when 
 /// progressing along the edge indicated by direction (S/E/N/W); 
 /// Left/lower vertex: S_base; Right/upper vertex: S_base + S_step
 static DenseMatrix<double> S_step; 
  
 /// Get opposite edge, e.g. Reflect_edge[N]=S 
 static Vector<int> Reflect_edge;

 /// \short Array of direction/quadrant adjacency scheme: 
 /// Is_adjacent(i_vertex_or_edge,j_quadrant): Is edge/vertex 
 /// adjacent to quadrant?
 static DenseMatrix<bool> Is_adjacent;

 /// \short Reflection scheme: Reflect(direction,quadrant): Get mirror 
 /// of quadrant in specified direction. E.g. Reflect(S,NE)=SE
 static  DenseMatrix<int> Reflect;

 /// \short Rotate coordinates: If North becomes NorthIs then direction 
 /// becomes Rotate(NorthIs,direction). E.g.  Rotate(E,NW)=NE;
 static DenseMatrix<int> Rotate;
 
 /// \short Angle betwen rotated coordinates: If old_direction becomes 
 /// new_direction then the angle between the axes (in anti-clockwise 
 /// direction is Rotate_angle(old_direction,new_direction); E.g. 
 /// Rotate_angle(E,N)=90;
 static DenseMatrix<int> Rotate_angle; 
 
 /// \short S_direct(direction,son_quadrant): The lower left corner 
 /// of son_quadrant has an offset of h/2 S_direct(direction,son_quadrant) 
 /// in the specified direction. E.g. S_direct(S,NE)=1 and  S_direct(S,NW)=0
 static DenseMatrix<int> S_direct;

};



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

 /// \short Vector giving the north equivalent of the neighbours:
 /// When viewed from the current quadtree's \c neighbour neighbour,
 /// our northern direction is the neighbour's North_equivalent[neighbour] 
 /// direction. If there's no rotation, this map contains the identify
 /// so that, e.g. \c North_equivalent[W]=N (read as: "in my Western
 /// neighbour, my North is its North"). If the western neighbour is rotated
 /// by 180 degrees relative to the current quadtree, say, we have
 /// \c North_equivalent[W]=S (read as: "in my Western
 /// neighbour, my North is its South"); etc.
 Vector<int> North_equivalent;

  public:

 ///Constructor for the (empty) root quadtree: Pass pointer to 
 /// associated object, a RefineableQElement<2>.
 QuadTreeRoot(RefineableElement* const &object_pt) : Tree(object_pt),
  QuadTree(object_pt), TreeRoot(object_pt)
  {
   
#ifdef PARANOID
   // Check that static member data has been setup
   if (!Static_data_has_been_setup)
    {
     std::string error_message =
      "Static member data hasn't been setup yet.\n";
     error_message +=
      "Call QuadTree::setup_static_data() before creating\n";
     error_message += "any QuadTreeRoots\n";
 
     throw OomphLibError(error_message,
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   using namespace QuadTreeNames;

   //Initialise the North equivalents of the neighbouring QuadTreeRoots
   North_equivalent.resize(27);

   North_equivalent[N] = N;
   North_equivalent[E] = N;
   North_equivalent[W] = N;
   North_equivalent[S] = N;
  }
 

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


 /// \short Return north equivalent of the neighbours in specified
 /// direction: When viewed from the current quadtree's \c neighbour neighbour,
 /// our northern direction is the neighbour's north_equivalent(neighbour) 
 /// direction. If there's no rotation, this map contains the identify
 /// so that, e.g. \c north_equivalent(W)=N (read as: "in my Western
 /// neighbour, my North is its North"). If the western neighbour is rotated
 /// by 180 degrees relative to the current quadtree, say, we have
 /// \c north_equivalent(W)=S (read as: "in my Western
 /// neighbour, my North is its South"); etc.
 int &north_equivalent(const int &neighbour)
  {
#ifdef PARANOID
   using namespace QuadTreeNames;   
   // Neighbour can only be located in N/S/E/W direction
   if((neighbour!=S)&&(neighbour!=N)&&(neighbour!=W)&&(neighbour!=E))
    {
     std::ostringstream error_message;
     error_message << "The neighbour can only be N,S,E,W, not" 
                   << Direct_string[neighbour] << std::endl;

     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return North_equivalent[neighbour];
  }


 /// \short If quadtree_root_pt is a neighbour, return the direction
 /// [N/S/E/W] in which it is found, otherwise return OMEGA
 int direction_of_neighbour(QuadTreeRoot* quadtree_root_pt)
  {
   using namespace QuadTreeNames;
   if(Neighbour_pt[N]==quadtree_root_pt) {return N;}
   if(Neighbour_pt[E]==quadtree_root_pt) {return E;}
   if(Neighbour_pt[S]==quadtree_root_pt) {return S;}
   if(Neighbour_pt[W]==quadtree_root_pt) {return W;}
   //If we get here, it's not a neighbour
   return OMEGA;
  }

};
 


//================================================================
/// A QuadTreeForest consists of a collection of QuadTreeRoots.
/// Each member tree can have neighbours to its S/W/N/E
/// and the orientation of their compasses can differ, allowing
/// for complex, unstructured meshes.
//=================================================================
class QuadTreeForest : public TreeForest
{ 
 public:

 /// Default constructor (empty and broken)
 QuadTreeForest()
  {
   //Throw an error
   throw OomphLibError(
    "Don't call an empty constructor for a QuadTreeForest object",
    OOMPH_CURRENT_FUNCTION,OOMPH_EXCEPTION_LOCATION);
  }
 
 /// \short Constructor: Pass vector of pointers to the roots of the
 /// constituent QuadTrees
 QuadTreeForest(Vector<TreeRoot* >& trees_pt);

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

 /// \short Destructor: Delete the constituent quadtrees (and thus 
 /// the objects associated with its non-leaf nodes!)
 virtual ~QuadTreeForest() {}
 
 /// \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 QuadTree::Max_neighbour_finding_tolerance.
 unsigned self_test();
    
private:

 /// Construct the rotation schemes
 void construct_north_equivalents();

 /// Construct the neighbour lookup scheme
 void find_neighbours();

 /// Return pointer to i-th root quadtree in this forest.
 /// (Performs a dynamic cast from the TreeRoot to a 
 /// QuadTreeRoot). 
 QuadTreeRoot* quadtree_pt(const unsigned &i)
  {return dynamic_cast<QuadTreeRoot*>(Trees_pt[i]);}
 
 /// \short Given the number i of the root quadtree in this forest, return
 /// pointer to its neighbour in the specified direction. NULL
 /// if neighbour doesn't exist. (This does the dynamic cast 
 /// from a TreeRoot to a QuadTreeRoot internally).
 QuadTreeRoot* quad_neigh_pt(const unsigned &i, const int &direction)
  {return dynamic_cast<QuadTreeRoot*>(Trees_pt[i]->neighbour_pt(direction));}

};

}

#endif