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-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//====================================================================
#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
  {
    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;
 
  //================================================================
  //=================================================================
  class OcTree : public virtual Tree
  {
  public:
    virtual ~OcTree() {}
 
 
    OcTree(const OcTree& dummy) = delete;
 
    void operator=(const OcTree&) = delete;
 
    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;
    }
 
    // lies between, i.e. the faces to which it is a common edge
    static Vector<int> faces_of_common_edge(const int& edge);
 
    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;
 
    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;
 
 
    unsigned self_test();
 
    static void setup_static_data();
 
 
    static void doc_face_neighbours(Vector<Tree*> forest_nodes_pt,
                                    std::ofstream& neighbours_file,
                                    std::ofstream& neighbours_txt_file,
                                    double& max_error);
 
 
    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);
 
    static int get_the_other_face(const unsigned& n1,
                                  const unsigned& n2,
                                  const unsigned& nnode1d,
                                  const int& face);
 
    static unsigned vertex_to_node_number(const int& vertex,
                                          const unsigned& nnode1d);
 
 
    static int node_number_to_vertex(const unsigned& n,
                                     const unsigned& nnode1d);
 
 
    static Vector<int> rotate(const int& new_up,
                              const int& new_right,
                              const Vector<int>& dir);
 
 
    static int rotate(const int& new_up, const int& new_right, const int& dir);
 
 
    static Vector<std::string> Direct_string;
 
    static Vector<int> Reflect_face;
 
    static Vector<int> Reflect_edge;
 
    static Vector<int> Reflect_vertex;
 
    static Vector<Vector<int>> Vertex_at_end_of_edge;
 
    static std::map<Vector<int>, int> Vector_to_direction;
 
    static Vector<Vector<int>> Direction_to_vector;
 
 
    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:
    OcTree()
    {
      throw OomphLibError("Don't call empty constructor for OcTree!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    OcTree(RefineableElement* const& object_pt) : Tree(object_pt) {}
 
 
    OcTree(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:
    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;
 
 
    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;
 
 
    bool edge_neighbour_is_face_neighbour(const int& edge,
                                          OcTree* edge_neighb_pt) const;
 
 
    static void construct_rotation_matrix(int& axis,
                                          int& angle,
                                          DenseMatrix<int>& mat);
 
    static void mult_mat_vect(const DenseMatrix<int>& mat,
                              const Vector<int>& vect1,
                              Vector<int>& vect2);
 
    static void mult_mat_mat(const DenseMatrix<int>& mat1,
                             const DenseMatrix<int>& mat2,
                             DenseMatrix<int>& mat3);
 
 
    static Vector<int> vertex_node_to_vector(const unsigned& n,
                                             const unsigned& nnode1d);
 
 
    static Vector<int> Cosi;
 
    static Vector<int> Sini;
 
 
    static DenseMatrix<bool> Is_adjacent;
 
    static DenseMatrix<int> Reflect;
 
    static DenseMatrix<int> Common_face;
 
    static Vector<std::string> Colour;
 
    static DenseMatrix<double> S_base;
 
    static DenseMatrix<double> S_steplo;
 
    static DenseMatrix<double> S_stephi;
 
    static DenseMatrix<double> S_directlo;
 
    static DenseMatrix<double> S_directhi;
 
    static DenseMatrix<double> S_base_edge;
 
    static DenseMatrix<double> S_step_edge;
 
    static DenseMatrix<double> S_direct_edge;
  };
 
 
  //===================================================================
  //==================================================================
  class OcTreeRoot : public virtual OcTree, public virtual TreeRoot
  {
  private:
    std::map<int, Vector<TreeRoot*>> Edge_neighbour_pt;
 
    std::map<TreeRoot*, int> Up_equivalent;
 
    std::map<TreeRoot*, int> Right_equivalent;
 
 
  public:
    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
    }
 
 
    OcTreeRoot(const OcTreeRoot& dummy) = delete;
 
    void operator=(const OcTreeRoot&) = delete;
 
    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];
    }
 
 
    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();
    }
 
    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);
      }
    }
 
 
    int up_equivalent(TreeRoot* tree_root_pt)
    {
      if (direction_of_neighbour(tree_root_pt) == OMEGA)
      {
        return OMEGA;
      }
      else
      {
        return Up_equivalent[tree_root_pt];
      }
    }
 
 
    void set_up_equivalent(TreeRoot* tree_root_pt, const int& dir)
    {
      Up_equivalent[tree_root_pt] = dir;
    }
 
 
    int right_equivalent(TreeRoot* tree_root_pt)
    {
      if (direction_of_neighbour(tree_root_pt) == OMEGA)
      {
        return OMEGA;
      }
      else
      {
        return Right_equivalent[tree_root_pt];
      }
    }
 
    void set_right_equivalent(TreeRoot* tree_root_pt, const int& dir)
    {
      Right_equivalent[tree_root_pt] = dir;
    }
 
    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;
    }
  };
 
 
 
 
  //================================================================
  //=================================================================
  class OcTreeForest : public TreeForest
  {
  public:
    OcTreeForest(Vector<TreeRoot*>& trees_pt);
 
    OcTreeForest()
    {
      throw OomphLibError("Don't call empty constructor for OcTreeForest!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    OcTreeForest(const OcTreeForest& dummy) = delete;
 
    void operator=(const OcTreeForest&) = delete;
 
    virtual ~OcTreeForest() {}
 
 
    void check_all_neighbours(DocInfo& doc_info);
 
    void open_hanging_node_files(DocInfo& doc_info,
                                 Vector<std::ofstream*>& output_stream);
 
    unsigned self_test();
 
    OcTreeRoot* octree_pt(const unsigned& i) const
    {
      return dynamic_cast<OcTreeRoot*>(Trees_pt[i]);
    }
 
    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));
    }
 
    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);
    }
 
    void construct_up_right_equivalents();
 
  private:
    void find_neighbours();
  };
 
} // namespace oomph
 
#endif