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// Copyright (C) 2006-2022 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 QuadTreeNames
  {
    enum
    {
      SW,
      SE,
      NW,
      NE,
      N,
      E,
      S,
      W,
      OMEGA = 26
    };
  }; // namespace QuadTreeNames
 
  // Forward class definition for class representing the root of a QuadTree
  class QuadTreeRoot;
 
  //================================================================
  //=================================================================
  class QuadTree : public virtual Tree
  {
  public:
    virtual ~QuadTree() {}
 
    QuadTree(const QuadTree& dummy) = delete;
 
    void operator=(const QuadTree&) = delete;
 
    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;
    }
 
    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;
 
    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;
 
    unsigned self_test();
 
    static void setup_static_data();
 
    static void doc_neighbours(Vector<Tree*> forest_nodes_pt,
                               std::ofstream& neighbours_file,
                               std::ofstream& neighbours_txt_file,
                               double& max_error);
 
 
    static Vector<std::string> Direct_string;
 
  protected:
    QuadTree()
    {
      throw OomphLibError(
        "Don't call an empty constructor for a QuadTree object",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    QuadTree(RefineableElement* const& object_pt) : Tree(object_pt) {}
 
    QuadTree(RefineableElement* const& object_pt,
             Tree* const& father_pt,
             const int& son_type)
      : Tree(object_pt, father_pt, son_type)
    {
    }
 
    static bool Static_data_has_been_setup;
 
 
  private:
    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;
 
    static Vector<std::string> Colour;
 
    static DenseMatrix<double> S_base;
 
    static DenseMatrix<double> S_step;
 
    static Vector<int> Reflect_edge;
 
    static DenseMatrix<bool> Is_adjacent;
 
    static DenseMatrix<int> Reflect;
 
    static DenseMatrix<int> Rotate;
 
    static DenseMatrix<int> Rotate_angle;
 
    static DenseMatrix<int> S_direct;
  };
 
 
  //===================================================================
  //==================================================================
  class QuadTreeRoot : public virtual QuadTree, public virtual TreeRoot
  {
  private:
    Vector<int> North_equivalent;
 
  public:
    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;
    }
 
 
    QuadTreeRoot(const QuadTreeRoot& dummy) = delete;
 
    void operator=(const QuadTreeRoot&) = delete;
 
 
    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];
    }
 
 
    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;
    }
  };
 
 
  //================================================================
  //=================================================================
  class QuadTreeForest : public TreeForest
  {
  public:
    QuadTreeForest()
    {
      // Throw an error
      throw OomphLibError(
        "Don't call an empty constructor for a QuadTreeForest object",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    QuadTreeForest(Vector<TreeRoot*>& trees_pt);
 
    QuadTreeForest(const QuadTreeForest& dummy) = delete;
 
    void operator=(const QuadTreeForest&) = delete;
 
    virtual ~QuadTreeForest() {}
 
    void check_all_neighbours(DocInfo& doc_info);
 
    void open_hanging_node_files(DocInfo& doc_info,
                                 Vector<std::ofstream*>& output_stream);
 
    unsigned self_test();
 
  private:
    void construct_north_equivalents();
 
    void find_neighbours();
 
    QuadTreeRoot* quadtree_pt(const unsigned& i)
    {
      return dynamic_cast<QuadTreeRoot*>(Trees_pt[i]);
    }
 
    QuadTreeRoot* quad_neigh_pt(const unsigned& i, const int& direction)
    {
      return dynamic_cast<QuadTreeRoot*>(Trees_pt[i]->neighbour_pt(direction));
    }
  };
 
} // namespace oomph
 
#endif