Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Static Protected Attributes | Private Member Functions | Static Private Member Functions | Static Private Attributes | List of all members
oomph::OcTree Class Reference

#include <octree.h>

+ Inheritance diagram for oomph::OcTree:

Public Member Functions

virtual ~OcTree ()
 
 OcTree (const OcTree &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const OcTree &)=delete
 Broken assignment operator. More...
 
Tree * construct_son (RefineableElement *const &object_pt, Tree *const &father_pt, const int &son_type)
 
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 ()
 
- Public Member Functions inherited from oomph::Tree
virtual ~Tree ()
 
 Tree (const Tree &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Tree &)=delete
 Broken assignment operator. More...
 
RefineableElement * object_pt () const
 
void flush_object ()
 Flush the object represented by the tree. More...
 
Tree * son_pt (const int &son_index) const
 
void set_son_pt (const Vector< Tree * > &son_pt)
 
unsigned nsons () const
 Return number of sons (zero if it's a leaf node) More...
 
void flush_sons ()
 Flush the sons. More...
 
TreeRoot *& root_pt ()
 Return pointer to root of the tree. More...
 
TreeRoot * root_pt () const
 Return pointer to root of the tree (const version) More...
 
template<class ELEMENT >
void split_if_required ()
 
template<class ELEMENT >
void p_refine_if_required (Mesh *&mesh_pt)
 
void merge_sons_if_required (Mesh *&mesh_pt)
 
void deactivate_object ()
 Call the RefineableElement's deactivate_element() function. More...
 
void traverse_all (Tree::VoidMemberFctPt member_function)
 
void traverse_all (Tree::VoidMeshPtArgumentMemberFctPt member_function, Mesh *&mesh_pt)
 
void traverse_all_but_leaves (Tree::VoidMemberFctPt member_function)
 
void traverse_leaves (Tree::VoidMemberFctPt member_function)
 
void traverse_leaves (Tree::VoidMeshPtArgumentMemberFctPt member_function, Mesh *&mesh_pt)
 
void stick_leaves_into_vector (Vector< Tree * > &)
 Traverse tree and stick pointers to leaf "nodes" (only) into Vector. More...
 
void stick_all_tree_nodes_into_vector (Vector< Tree * > &)
 Traverse and stick pointers to all "nodes" into Vector. More...
 
int son_type () const
 Return son type. More...
 
bool is_leaf ()
 Return true if the tree is a leaf node. More...
 
Tree * father_pt () const
 Return pointer to father: NULL if it's a root node. More...
 
void set_father_pt (Tree *const &father_pt)
 Set the father. More...
 
unsigned level () const
 Return the level of the Tree (root=0) More...
 

Static Public Member Functions

static Vector< int > faces_of_common_edge (const int &edge)
 Function that, given an edge, returns the two faces on which it. More...
 
static void setup_static_data ()
 Setup the static data, rotation and reflection schemes, etc. More...
 
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 Public Member Functions inherited from oomph::Tree
static double & max_neighbour_finding_tolerance ()
 

Static Public Attributes

static Vector< std::string > Direct_string
 Translate (enumerated) directions into strings. More...
 
static Vector< int > Reflect_face
 Get opposite face, e.g. Reflect_face[L]=R. More...
 
static Vector< int > Reflect_edge
 Get opposite edge, e.g. Reflect_edge[DB]=UF. More...
 
static Vector< int > Reflect_vertex
 Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF. More...
 
static Vector< Vector< int > > Vertex_at_end_of_edge
 Map of vectors containing the two vertices for each edge. More...
 
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
 
- Static Public Attributes inherited from oomph::Tree
static const int OMEGA = 26
 Default value for an unassigned neighbour. More...
 

Protected Member Functions

 OcTree ()
 Default constructor (empty and broken) More...
 
 OcTree (RefineableElement *const &object_pt)
 
 OcTree (RefineableElement *const &object_pt, Tree *const &father_pt, const int &son_type)
 
- Protected Member Functions inherited from oomph::Tree
 Tree ()
 Default constructor (empty and broken) More...
 
 Tree (RefineableElement *const &object_pt)
 
 Tree (RefineableElement *const &object_pt, Tree *const &father_pt, const int &son_type)
 

Static Protected Attributes

static bool Static_data_has_been_setup = false
 Bool indicating that static member data has been setup. More...
 
- Static Protected Attributes inherited from oomph::Tree
static double Max_neighbour_finding_tolerance = 1.0e-14
 

Private Member Functions

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 Private Member Functions

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)
 Helper function: Performs the operation : vect2 = mat*vect1. More...
 
static void mult_mat_mat (const DenseMatrix< int > &mat1, const DenseMatrix< int > &mat2, DenseMatrix< int > &mat3)
 Helper function: Performs the operation : mat3=mat1*mat2. More...
 
static Vector< int > vertex_node_to_vector (const unsigned &n, const unsigned &nnode1d)
 

Static Private Attributes

static Vector< int > Cosi
 Entry in rotation matrix: cos(i*90) More...
 
static Vector< int > Sini
 Entry in rotation matrix sin(i*90) More...
 
static DenseMatrix< bool > Is_adjacent
 
static DenseMatrix< int > Reflect
 
static DenseMatrix< int > Common_face
 
static Vector< std::string > Colour
 Colours for neighbours in various directions. More...
 
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
 

Additional Inherited Members

- Public Types inherited from oomph::Tree
typedef void(Tree::* VoidMemberFctPt) ()
 Function pointer to argument-free void Tree member function. More...
 
typedef void(Tree::* VoidMeshPtArgumentMemberFctPt) (Mesh *&mesh_pt)
 
- Protected Attributes inherited from oomph::Tree
TreeRoot * Root_pt
 Pointer to the root of the tree. More...
 
Tree * Father_pt
 Pointer to the Father of the Tree. More...
 
Vector< Tree * > Son_pt
 Vector of pointers to the sons of the Tree. More...
 
int Level
 Level of the Tree (level 0 = root) More...
 
int Son_type
 Son type (e.g. SW/SE/NW/NE in a quadtree) More...
 
RefineableElement * Object_pt
 Pointer to the object represented by the tree. More...
 

Detailed Description

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 \( {\bf s} \in [-1,1]^3 \).

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.

Constructor & Destructor Documentation

◆ ~OcTree()

virtual oomph::OcTree::~OcTree ( )
inlinevirtual

Destructor. Note: Deleting a octree also deletes the objects associated with all non-leaf nodes!

118 {}

◆ OcTree() [1/4]

oomph::OcTree::OcTree ( const OcTree &  dummy)
delete

Broken copy constructor.

◆ OcTree() [2/4]

oomph::OcTree::OcTree ( )
inlineprotected

Default constructor (empty and broken)

383  {
384  throw OomphLibError("Don't call empty constructor for OcTree!",
385  OOMPH_CURRENT_FUNCTION,
386  OOMPH_EXCEPTION_LOCATION);
387  }
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by construct_son().

◆ OcTree() [3/4]

oomph::OcTree::OcTree ( RefineableElement *const &  object_pt)
inlineprotected

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.

395 : Tree(object_pt) {}
oomph::Tree::object_pt
RefineableElement * object_pt() const
Definition: tree.h:88
oomph::Tree::Tree
Tree()
Default constructor (empty and broken)
Definition: tree.h:266

◆ OcTree() [4/4]

oomph::OcTree::OcTree ( RefineableElement *const &  object_pt,
Tree *const &  father_pt,
const int &  son_type 
)
inlineprotected

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.

407  : Tree(object_pt, father_pt, son_type)
408  {
409  }
oomph::Tree::father_pt
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:235
oomph::Tree::son_type
int son_type() const
Return son type.
Definition: tree.h:214

Member Function Documentation

◆ construct_rotation_matrix()

void oomph::OcTree::construct_rotation_matrix ( int &  axis,
int &  angle,
DenseMatrix< int > &  mat 
)
staticprivate

This constructs the rotation matrix of the rotation around the axis axis with an angle of angle*90

Build the rotation matrix for a rotation around the axis axis of an angle angle*90

612  {
613  using namespace OcTreeNames;
614  int a = 0, b = 0, c = 0, i, j;
615 
616  switch (axis)
617  {
618  case R:
619  a = 1;
620  b = 2;
621  c = 0;
622  break;
623  case U:
624  a = 2;
625  b = 0;
626  c = 1;
627  break;
628  case F:
629  a = 0;
630  b = 1;
631  c = 2;
632  break;
633  default:
634  // Object to create error message
635  std::ostringstream error_message_stream;
636 
637  // Create the error message
638  error_message_stream << "Bad axis (" << axis << "). Expected "
639  << OcTreeNames::R << ", " << OcTreeNames::U
640  << " or " << OcTreeNames::F << "." << std::endl;
641 
642  // Throw the error message
643  throw OomphLibError(error_message_stream.str(),
644  OOMPH_CURRENT_FUNCTION,
645  OOMPH_EXCEPTION_LOCATION);
646  }
647  for (i = 0; i < 3; i++)
648  {
649  for (j = 0; j < 3; j++)
650  {
651  mat(i, j) = 0;
652  }
653  }
654  mat(a, a) = Cosi[angle];
655  mat(b, b) = Cosi[angle];
656  mat(a, b) = -Sini[angle];
657  mat(b, a) = Sini[angle];
658  mat(c, c) = 1;
659  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
mat
Eigen::SparseMatrix< double > mat
Definition: EigenUnitTest.cpp:10
R
@ R
Definition: StatisticsVector.h:21
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::OcTree::Sini
static Vector< int > Sini
Entry in rotation matrix sin(i*90)
Definition: octree.h:523
oomph::OcTree::Cosi
static Vector< int > Cosi
Entry in rotation matrix: cos(i*90)
Definition: octree.h:520
a
const Scalar * a
Definition: level2_cplx_impl.h:32
Jeffery_Solution::angle
double angle(const double &t)
Angular position as a function of time t.
Definition: jeffery_orbit.cc:98
RachelsAdvectionDiffusion::U
double U
Swimming speed.
Definition: two_d_variable_diff_adapt.cc:53
calibrate.c
int c
Definition: calibrate.py:100
oomph::OcTreeNames::R
@ R
Definition: octree.h:70
oomph::OcTreeNames::F
@ F
Definition: octree.h:74
oomph::OcTreeNames::U
@ U
Definition: octree.h:72
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References a, Jeffery_Solution::angle(), b, calibrate::c, Cosi, oomph::OcTreeNames::F, i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, oomph::OcTreeNames::R, Sini, RachelsAdvectionDiffusion::U, and oomph::OcTreeNames::U.

Referenced by rotate().

◆ construct_son()

Tree* oomph::OcTree::construct_son ( RefineableElement *const &  object_pt,
Tree *const &  father_pt,
const int &  son_type 
)
inlinevirtual

Overload the function construct_son to ensure that the son is a specific OcTree and not a general Tree.

Implements oomph::Tree.

132  {
133  OcTree* temp_oc_pt = new OcTree(object_pt, father_pt, son_type);
134  return temp_oc_pt;
135  }
oomph::OcTree::OcTree
OcTree()
Default constructor (empty and broken)
Definition: octree.h:382

References oomph::Tree::father_pt(), oomph::Tree::object_pt(), OcTree(), and oomph::Tree::son_type().

◆ doc_face_neighbours()

void oomph::OcTree::doc_face_neighbours ( Vector< Tree * >  forest_nodes_pt,
std::ofstream &  neighbours_file,
std::ofstream &  neighbours_txt_file,
double &  max_error 
)
static

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.

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 neighbouring element. If the two filestreams are closed, output is suppressed. (Static function.)

4279  {
4280  using namespace OcTreeNames;
4281 
4282  int diff_level;
4283  int face = OMEGA;
4284  bool in_neighbouring_tree;
4285 
4286  Vector<double> s(3);
4287  Vector<double> x(3);
4288 
4289  Vector<double> s_sw(3);
4290  Vector<double> s_ne(3);
4291  Vector<unsigned> translate_s(3);
4292 
4293  Vector<double> x_small(3);
4294  Vector<double> x_large(3);
4295 
4296 
4297  // Initialise error in vertex positions
4298  max_error = 0.0;
4299 
4300  // Loop over all elements
4301  // ----------------------
4302  unsigned long num_nodes = forest_nodes_pt.size();
4303 
4304  for (unsigned long i = 0; i < num_nodes; i++)
4305  {
4306  // Doc the element itself
4307  OcTree* el_pt = dynamic_cast<OcTree*>(forest_nodes_pt[i]);
4308 
4309  // If the object is incomplete omit
4310  if (el_pt->object_pt()->nodes_built())
4311  {
4312  // Print it
4313  if (neighbours_file.is_open())
4314  {
4315  neighbours_file << "#---------------------------------" << std::endl;
4316  neighbours_file << "#The element itself: " << i << std::endl;
4317  neighbours_file << "#---------------------------------" << std::endl;
4318  dynamic_cast<RefineableQElement<3>*>(el_pt->object_pt())
4319  ->output_corners(neighbours_file, "BLACK");
4320  }
4321 
4322  // Loop over directions to find neighbours
4323  // ----------------------------------------
4324  for (int direction = L; direction <= F; direction++)
4325  {
4326  // Initialise difference in levels and coordinate offset
4327  diff_level = 0;
4328 
4329  // Find greater-or-equal-sized neighbour...
4330  OcTree* neighb_pt = el_pt->gteq_face_neighbour(direction,
4331  translate_s,
4332  s_sw,
4333  s_ne,
4334  face,
4335  diff_level,
4336  in_neighbouring_tree);
4337 
4338  // If neighbour exists and nodes are created
4339  if ((neighb_pt != 0) && (neighb_pt->object_pt()->nodes_built()))
4340  {
4341  // Doc neighbour stats
4342  if (neighbours_txt_file.is_open())
4343  {
4344  neighbours_txt_file
4345  << Direct_string[direction] << " neighbour of "
4346  << el_pt->object_pt()->number() << " is "
4347  << neighb_pt->object_pt()->number() << " diff_level "
4348  << diff_level << ". Inside the neighbour the face is "
4349  << Direct_string[face] << std::endl;
4350  }
4351 
4352  // Plot neighbour in the appropriate colour
4353  if (neighbours_file.is_open())
4354  {
4355  neighbours_file << "#---------------------------------"
4356  << std::endl;
4357  neighbours_file
4358  << "#Neighbour element: " << Direct_string[direction]
4359  << "\n#---------------------------------" << std::endl;
4360  dynamic_cast<RefineableQElement<3>*>(neighb_pt->object_pt())
4361  ->output_corners(neighbours_file, Colour[direction]);
4362  }
4363 
4364  // Check that local coordinates in the larger element
4365  // lead to the same spatial point as the node vertices
4366  // in the current element
4367  if (neighbours_file.is_open())
4368  {
4369  neighbours_file << "ZONE I=2 C=" << Colour[direction]
4370  << std::endl;
4371  }
4372 
4373  // "South west" vertex in the interfacial face
4374  //--------------------------------------------
4375 
4376  // Get coordinates in large (neighbour) element
4377  s[0] = s_sw[0];
4378  s[1] = s_sw[1];
4379  s[2] = s_sw[2];
4380  neighb_pt->object_pt()->get_x(s, x_large);
4381 
4382  // Get coordinates in small element
4383  Vector<double> s(3);
4384  s[0] = S_base(0, direction);
4385  s[1] = S_base(1, direction);
4386  s[2] = S_base(2, direction);
4387  el_pt->object_pt()->get_x(s, x_small);
4388 
4389  // Need to exclude periodic nodes from this check
4390  // There can only be periodic nodes if we have moved into the
4391  // neighbour
4392  bool is_periodic = false;
4393  if (in_neighbouring_tree)
4394  {
4395  // is the node periodic
4396  is_periodic = el_pt->root_pt()->is_neighbour_periodic(direction);
4397  }
4398 
4399  double error = 0.0;
4400  // Only bother to calculate the error if the node is NOT periodic
4401  if (is_periodic == false)
4402  {
4403  error = sqrt(pow(x_small[0] - x_large[0], 2) +
4404  pow(x_small[1] - x_large[1], 2) +
4405  pow(x_small[2] - x_large[2], 2));
4406  }
4407 
4408  if (neighbours_txt_file.is_open())
4409  {
4410  neighbours_txt_file << "Error (1) " << error << std::endl;
4411  }
4412 
4413  if (std::fabs(error) > max_error)
4414  {
4415  max_error = std::fabs(error);
4416  }
4417 
4418  // Check error and doc mismatch if required
4419  bool stop = false;
4420  std::ofstream mismatch_file;
4421  if (std::fabs(error) > Max_neighbour_finding_tolerance)
4422  {
4423  stop = true;
4424  mismatch_file.open("mismatch.dat");
4425  mismatch_file << "ZONE" << std::endl;
4426  mismatch_file << x_large[0] << " " << x_large[1] << " "
4427  << x_large[2] << " 2 \n";
4428  mismatch_file << x_small[0] << " " << x_small[1] << " "
4429  << x_small[2] << " 3 \n";
4430  }
4431 
4432  if (neighbours_file.is_open())
4433  {
4434  neighbours_file << "#SOUTH WEST: " << std::endl;
4435  neighbours_file << x_large[0] << " " << x_large[1] << " "
4436  << x_large[2] << " 40 \n";
4437  }
4438 
4439 
4440  // "North east" vertex in the interfacial face
4441  //--------------------------------------------
4442 
4443  // Get coordinates in large (neighbour) element
4444  s[0] = s_ne[0];
4445  s[1] = s_ne[1];
4446  s[2] = s_ne[2];
4447  neighb_pt->object_pt()->get_x(s, x_large);
4448 
4449  // Get coordinates in small element
4450  s[0] = S_base(0, direction) + S_steplo(0, direction) +
4451  S_stephi(0, direction);
4452  s[1] = S_base(1, direction) + S_steplo(1, direction) +
4453  S_stephi(1, direction);
4454  s[2] = S_base(2, direction) + S_steplo(2, direction) +
4455  S_stephi(2, direction);
4456  el_pt->object_pt()->get_x(s, x_small);
4457 
4458  error = 0.0;
4459  // Only do this if we are NOT periodic
4460  if (is_periodic == false)
4461  {
4462  error = sqrt(pow(x_small[0] - x_large[0], 2) +
4463  pow(x_small[1] - x_large[1], 2) +
4464  pow(x_small[2] - x_large[2], 2));
4465  }
4466 
4467  // Output
4468  if (neighbours_file.is_open())
4469  {
4470  neighbours_file << "#NORTH EAST: " << std::endl;
4471  neighbours_file << x_large[0] << " " << x_large[1] << " "
4472  << x_large[2] << " 80 \n";
4473  }
4474 
4475  if (neighbours_txt_file.is_open())
4476  {
4477  neighbours_txt_file << "Error (2) " << error << std::endl;
4478  }
4479 
4480  if (std::fabs(error) > max_error)
4481  {
4482  max_error = std::fabs(error);
4483  }
4484 
4485  // Check error and doc mismatch if required
4486  if (std::fabs(error) > Max_neighbour_finding_tolerance)
4487  {
4488  stop = true;
4489  }
4490  if (stop)
4491  {
4492  if (!mismatch_file.is_open())
4493  {
4494  mismatch_file.open("mismatch.dat");
4495  }
4496  mismatch_file << "ZONE" << std::endl;
4497  mismatch_file << x_large[0] << " " << x_large[1] << " "
4498  << x_large[2] << " 2 " << std::fabs(error)
4499  << " \n";
4500  mismatch_file << x_small[0] << " " << x_small[1] << " "
4501  << x_small[2] << " 3 " << std::fabs(error)
4502  << " \n";
4503  mismatch_file.close();
4504  pause("Error");
4505  }
4506  }
4507 
4508  // If neighbour does not exist: Insert blank zones into file
4509  // so that tecplot can find six neighbours for every element
4510  else
4511  {
4512  if (neighbours_file.is_open())
4513  {
4514  neighbours_file
4515  << "#---------------------------------\n"
4516  << "# No neighbour in direction: " << Direct_string[direction]
4517  << "\n"
4518  << "#---------------------------------" << std::endl;
4519 
4520  dynamic_cast<RefineableQElement<3>*>(el_pt->object_pt())
4521  ->output_corners(neighbours_file, "WHITE");
4522  neighbours_file << "ZONE I=2 \n";
4523  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4524  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4525  }
4526  }
4527 
4528  if (neighbours_file.is_open())
4529  {
4530  neighbours_file << std::endl << std::endl << std::endl;
4531  }
4532  }
4533  } // End of case when element can be documented
4534  }
4535  }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
L
MatrixXd L
Definition: LLT_example.cpp:6
oomph::OcTree::Colour
static Vector< std::string > Colour
Colours for neighbours in various directions.
Definition: octree.h:540
oomph::OcTree::S_base
static DenseMatrix< double > S_base
Definition: octree.h:544
oomph::OcTree::Direct_string
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:329
oomph::OcTree::S_stephi
static DenseMatrix< double > S_stephi
Definition: octree.h:565
oomph::OcTree::S_steplo
static DenseMatrix< double > S_steplo
Definition: octree.h:558
oomph::Tree::OMEGA
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
oomph::Tree::Max_neighbour_finding_tolerance
static double Max_neighbour_finding_tolerance
Definition: tree.h:313
s
RealScalar s
Definition: level1_cplx_impl.h:130
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
MeshRefinement::max_error
double max_error
Definition: MortaringCantileverCompareToNonMortaring.cpp:188
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
calibrate.error
int error
Definition: calibrate.py:297
oomph::pause
void pause(std::string message)
Pause and display message.
Definition: oomph_utilities.cc:1265
plotDoE.x
list x
Definition: plotDoE.py:28

References Colour, Direct_string, calibrate::error, oomph::OcTreeNames::F, boost::multiprecision::fabs(), oomph::FiniteElement::get_x(), gteq_face_neighbour(), i, oomph::TreeRoot::is_neighbour_periodic(), L, MeshRefinement::max_error, oomph::Tree::Max_neighbour_finding_tolerance, oomph::RefineableElement::nodes_built(), oomph::RefineableElement::number(), oomph::Tree::object_pt(), oomph::Tree::OMEGA, oomph::pause(), Eigen::bfloat16_impl::pow(), oomph::Tree::root_pt(), s, S_base, S_stephi, S_steplo, sqrt(), and plotDoE::x.

Referenced by oomph::OcTreeForest::check_all_neighbours(), self_test(), and oomph::OcTreeForest::self_test().

◆ doc_true_edge_neighbours()

void oomph::OcTree::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

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.

/////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// 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 neighbouring element. If the two filestreams are closed, output is suppressed. (Static function).

4557  {
4558  using namespace OcTreeNames;
4559 
4560  int diff_level;
4561  int edge = OMEGA;
4562 
4563  Vector<double> s(3);
4564  Vector<double> x(3);
4565 
4566  Vector<double> s_lo(3);
4567  Vector<double> s_hi(3);
4568  Vector<unsigned> translate_s(3);
4569 
4570  Vector<double> x_small(3);
4571  Vector<double> x_large(3);
4572 
4573 
4574  // Initialise error in vertex positions
4575  max_error = 0.0;
4576 
4577  // Loop over all elements
4578  // ----------------------
4579  unsigned long num_nodes = forest_nodes_pt.size();
4580 
4581  for (unsigned long i = 0; i < num_nodes; i++)
4582  {
4583  // Doc the element itself
4584  OcTree* el_pt = dynamic_cast<OcTree*>(forest_nodes_pt[i]);
4585 
4586  // If the object is incomplete omit
4587  if (el_pt->object_pt()->nodes_built())
4588  {
4589  // Print it
4590  if (neighbours_file.is_open())
4591  {
4592  neighbours_file << "#---------------------------------" << std::endl;
4593  neighbours_file << "# The element itself: " << i << std::endl;
4594  neighbours_file << "#---------------------------------" << std::endl;
4595  dynamic_cast<RefineableQElement<3>*>(el_pt->object_pt())
4596  ->output_corners(neighbours_file, "BLACK");
4597  }
4598 
4599  // Loop over directions to find edge neighbours
4600  // --------------------------------------------
4601  for (int direction = LB; direction <= UF; direction++)
4602  {
4603  // Initialise difference in levels
4604  diff_level = 0;
4605 
4606  // For now simply doc the zero-th edge neighbour (if any)
4607  unsigned i_root_edge_neighbour = 0;
4608  unsigned nroot_edge_neighbour = 0;
4609 
4610  // Find greater-or-equal-sized edge neighbour...
4611  OcTree* neighb_pt =
4612  el_pt->gteq_true_edge_neighbour(direction,
4613  i_root_edge_neighbour,
4614  nroot_edge_neighbour,
4615  translate_s,
4616  s_lo,
4617  s_hi,
4618  edge,
4619  diff_level);
4620 
4621  // If neighbour exists and nodes are created
4622  if ((neighb_pt != 0) && (neighb_pt->object_pt()->nodes_built()))
4623  {
4624  // Doc neighbour stats
4625  if (neighbours_txt_file.is_open())
4626  {
4627  neighbours_txt_file
4628  << Direct_string[direction] << " neighbour of "
4629  << el_pt->object_pt()->number() << " is "
4630  << neighb_pt->object_pt()->number() << " diff_level "
4631  << diff_level << ". Inside the neighbour the edge is "
4632  << Direct_string[edge] << std::endl;
4633  }
4634 
4635  // Plot neighbour in the appropriate colour
4636  if (neighbours_file.is_open())
4637  {
4638  neighbours_file
4639  << "#---------------------------------"
4640  << "\n# Neighbour element: " << Direct_string[direction]
4641  << "\n#---------------------------------" << std::endl;
4642  dynamic_cast<RefineableQElement<3>*>(neighb_pt->object_pt())
4643  ->output_corners(neighbours_file, Colour[direction]);
4644  }
4645 
4646  // Check that local coordinates in the larger element
4647  // lead to the same spatial point as the node vertices
4648  // in the current element
4649  if (neighbours_file.is_open())
4650  {
4651  neighbours_file << "ZONE I=2 C=" << Colour[direction]
4652  << std::endl;
4653  }
4654 
4655  // "Low" vertex in the edge
4656  //-------------------------
4657  // Get coordinates in large (neighbour) element
4658  s[0] = s_lo[0];
4659  s[1] = s_lo[1];
4660  s[2] = s_lo[2];
4661  neighb_pt->object_pt()->get_x(s, x_large);
4662 
4663  // Get coordinates in small element
4664  Vector<double> s(3);
4665  s[0] = S_base_edge(0, direction);
4666  s[1] = S_base_edge(1, direction);
4667  s[2] = S_base_edge(2, direction);
4668  el_pt->object_pt()->get_x(s, x_small);
4669 
4670  // Need to exclude periodic nodes from this check
4671  // There can only be periodic nodes if we have moved into the
4672  // neighbour
4673  bool is_periodic = false;
4674 
4675  // Get the faces on which the edge lies
4676  Vector<int> faces_attached_to_edge =
4677  faces_of_common_edge(direction);
4678 
4679  // Get the number of entries in the vector
4680  unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
4681 
4682  // Loop over the faces
4683  for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
4684  i_face++)
4685  {
4686  // Is the node periodic in the face direction?
4687  is_periodic = el_pt->root_pt()->is_neighbour_periodic(
4688  faces_attached_to_edge[i_face]);
4689 
4690  // Check if the edge is periodic in the i_face-th face direction
4691  if (is_periodic)
4692  {
4693  // We're done!
4694  break;
4695  }
4696  } // for (unsigned
4697  // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
4698 
4699  double error = 0.0;
4700  // Only bother to calculate the error if the node is NOT periodic
4701  if (is_periodic == false)
4702  {
4703  error = sqrt(pow(x_small[0] - x_large[0], 2) +
4704  pow(x_small[1] - x_large[1], 2) +
4705  pow(x_small[2] - x_large[2], 2));
4706  }
4707 
4708  if (std::fabs(error) > max_error)
4709  {
4710  max_error = std::fabs(error);
4711  }
4712 
4713  if (neighbours_txt_file.is_open())
4714  {
4715  neighbours_txt_file << "Error (1) " << error << std::endl;
4716  }
4717 
4718  // Check error and doc mismatch if required
4719  bool stop = false;
4720  std::ofstream mismatch_file;
4721  if (std::fabs(error) > Max_neighbour_finding_tolerance)
4722  {
4723  stop = true;
4724  mismatch_file.open("mismatch.dat");
4725  mismatch_file << "ZONE" << std::endl;
4726  mismatch_file << x_large[0] << " " << x_large[1] << " "
4727  << x_large[2] << " 2 \n";
4728  mismatch_file << x_small[0] << " " << x_small[1] << " "
4729  << x_small[2] << " 3 \n";
4730  }
4731 
4732  if (neighbours_file.is_open())
4733  {
4734  neighbours_file << "# LOW VERTEX: " << std::endl;
4735  neighbours_file << x_large[0] << " " << x_large[1] << " "
4736  << x_large[2] << " 40 \n";
4737  }
4738 
4739 
4740  // "High" vertex in the edge
4741  //--------------------------
4742  // Get coordinates in large (neighbour) element
4743  s[0] = s_hi[0];
4744  s[1] = s_hi[1];
4745  s[2] = s_hi[2];
4746  neighb_pt->object_pt()->get_x(s, x_large);
4747 
4748  // Get coordinates in small element
4749  s[0] = S_base_edge(0, direction) + S_step_edge(0, direction);
4750  s[1] = S_base_edge(1, direction) + S_step_edge(1, direction);
4751  s[2] = S_base_edge(2, direction) + S_step_edge(2, direction);
4752  el_pt->object_pt()->get_x(s, x_small);
4753 
4754  // Output
4755  if (neighbours_file.is_open())
4756  {
4757  neighbours_file << "# HI VERTEX: " << std::endl;
4758  neighbours_file << x_large[0] << " " << x_large[1] << " "
4759  << x_large[2] << " 80 \n";
4760  }
4761 
4762  // Reset the error value
4763  error = 0.0;
4764 
4765  // Only do this if we are NOT periodic
4766  if (is_periodic == false)
4767  {
4768  error = sqrt(pow(x_small[0] - x_large[0], 2) +
4769  pow(x_small[1] - x_large[1], 2) +
4770  pow(x_small[2] - x_large[2], 2));
4771  }
4772 
4773  if (neighbours_txt_file.is_open())
4774  {
4775  neighbours_txt_file << "Error (2) " << error << std::endl;
4776  }
4777 
4778  if (std::fabs(error) > max_error)
4779  {
4780  max_error = std::fabs(error);
4781  }
4782 
4783  // Check error and doc mismatch if required
4784  if (std::fabs(error) > Max_neighbour_finding_tolerance)
4785  {
4786  stop = true;
4787  }
4788  if (stop)
4789  {
4790  if (!mismatch_file.is_open())
4791  {
4792  mismatch_file.open("mismatch.dat");
4793  }
4794  mismatch_file << "ZONE" << std::endl;
4795  mismatch_file << x_large[0] << " " << x_large[1] << " "
4796  << x_large[2] << " 2 \n";
4797  mismatch_file << x_small[0] << " " << x_small[1] << " "
4798  << x_small[2] << " 3 \n";
4799  mismatch_file.close();
4800  pause("Error");
4801  }
4802  }
4803 
4804  // If neighbour does not exist: Insert blank zones into file
4805  // so that tecplot can find twelve neighbours for every element
4806  else
4807  {
4808  if (neighbours_file.is_open())
4809  {
4810  neighbours_file << "#---------------------------------"
4811  << std::endl;
4812  neighbours_file
4813  << "# No neighbour in direction: " << Direct_string[direction]
4814  << std::endl;
4815  neighbours_file << "#---------------------------------"
4816  << std::endl;
4817 
4818  dynamic_cast<RefineableQElement<3>*>(el_pt->object_pt())
4819  ->output_corners(neighbours_file, "WHITE");
4820  neighbours_file << "ZONE I=2 \n";
4821  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4822  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4823 
4824 
4825  // Doc edge for which no neighbour exists but only
4826  // for the smallest elements
4827  if (el_pt->is_leaf())
4828  {
4829  // Get start coordinates in small element
4830  Vector<double> s(3), x_start(3), x_end(3);
4831  s[0] = S_base_edge(0, direction);
4832  s[1] = S_base_edge(1, direction);
4833  s[2] = S_base_edge(2, direction);
4834  el_pt->object_pt()->get_x(s, x_start);
4835 
4836 
4837  // Get coordinates in small element
4838  s[0] = S_base_edge(0, direction) + S_step_edge(0, direction);
4839  s[1] = S_base_edge(1, direction) + S_step_edge(1, direction);
4840  s[2] = S_base_edge(2, direction) + S_step_edge(2, direction);
4841  el_pt->object_pt()->get_x(s, x_end);
4842 
4843  no_true_edge_file << "ZONE I=2" << std::endl;
4844  no_true_edge_file << x_start[0] << " " << x_start[1] << " "
4845  << x_start[2] << " " << std::endl;
4846  no_true_edge_file << x_end[0] << " " << x_end[1] << " "
4847  << x_end[2] << " " << std::endl;
4848  }
4849  }
4850 
4851  if (neighbours_file.is_open())
4852  {
4853  neighbours_file << std::endl << std::endl << std::endl;
4854  }
4855  }
4856  } // End of case when element can be documented
4857  }
4858  }
4859  }
oomph::OcTree::S_step_edge
static DenseMatrix< double > S_step_edge
Definition: octree.h:593
oomph::OcTree::faces_of_common_edge
static Vector< int > faces_of_common_edge(const int &edge)
Function that, given an edge, returns the two faces on which it.
Definition: octree.cc:268
oomph::OcTree::S_base_edge
static DenseMatrix< double > S_base_edge
Definition: octree.h:585
oomph::OcTreeNames::UF
@ UF
Definition: octree.h:68
oomph::OcTreeNames::LB
@ LB
Definition: octree.h:57

References Colour, Direct_string, calibrate::error, boost::multiprecision::fabs(), faces_of_common_edge(), oomph::FiniteElement::get_x(), gteq_true_edge_neighbour(), i, oomph::Tree::is_leaf(), oomph::TreeRoot::is_neighbour_periodic(), oomph::OcTreeNames::LB, MeshRefinement::max_error, oomph::Tree::Max_neighbour_finding_tolerance, oomph::RefineableElement::nodes_built(), oomph::RefineableElement::number(), oomph::Tree::object_pt(), oomph::Tree::OMEGA, oomph::pause(), Eigen::bfloat16_impl::pow(), oomph::Tree::root_pt(), s, S_base_edge, S_step_edge, sqrt(), oomph::OcTreeNames::UF, and plotDoE::x.

Referenced by oomph::OcTreeForest::check_all_neighbours(), self_test(), and oomph::OcTreeForest::self_test().

◆ edge_neighbour_is_face_neighbour()

bool oomph::OcTree::edge_neighbour_is_face_neighbour ( const int &  edge,
OcTree *  edge_neigh_pt 
) const
private

Is the edge neighbour (for edge "edge") specified via the pointer also a face neighbour for one of the two adjacent faces?

2771  {
2772 #ifdef PARANOID
2773  // No paranoid check needed -- the default for the switch statement
2774  // catches illegal values for edge
2775 #endif
2776 
2777 
2778  // Catch stupid case: Null doesn't have a face neighbour...
2779  if (edge_neigh_pt == 0)
2780  {
2781  return false;
2782  }
2783 
2784  using namespace OcTreeNames;
2785 
2786  // Auxiliary variables
2787  int face;
2788  Vector<unsigned> translate_s(3);
2789  Vector<double> s_sw(3);
2790  Vector<double> s_ne(3);
2791  int reflected_face;
2792  int diff_level;
2793  bool in_neighbouring_tree = false;
2794 
2795  OcTree* face_neigh_pt = 0;
2796 
2797  switch (edge)
2798  {
2799  case LB:
2800 
2801  // Get first face neighbour
2802  face = L;
2803  face_neigh_pt = gteq_face_neighbour(face,
2804  translate_s,
2805  s_sw,
2806  s_ne,
2807  reflected_face,
2808  diff_level,
2809  in_neighbouring_tree);
2810 
2811  // Check if they agree...
2812  if (face_neigh_pt != 0)
2813  {
2814  if (face_neigh_pt == edge_neigh_pt)
2815  {
2816  return true;
2817  }
2818  }
2819 
2820  // Get second face neighbour
2821  face = B;
2822  face_neigh_pt = gteq_face_neighbour(face,
2823  translate_s,
2824  s_sw,
2825  s_ne,
2826  reflected_face,
2827  diff_level,
2828  in_neighbouring_tree);
2829 
2830  // Check if they agree...
2831  if (face_neigh_pt != 0)
2832  {
2833  if (face_neigh_pt == edge_neigh_pt)
2834  {
2835  return true;
2836  }
2837  }
2838 
2839  break;
2840 
2841 
2842  case RB:
2843 
2844 
2845  // Get first face neighbour
2846  face = R;
2847  face_neigh_pt = gteq_face_neighbour(face,
2848  translate_s,
2849  s_sw,
2850  s_ne,
2851  reflected_face,
2852  diff_level,
2853  in_neighbouring_tree);
2854 
2855  // Check if they agree...
2856  if (face_neigh_pt != 0)
2857  {
2858  if (face_neigh_pt == edge_neigh_pt)
2859  {
2860  return true;
2861  }
2862  }
2863 
2864  // Get second face neighbour
2865  face = B;
2866  face_neigh_pt = gteq_face_neighbour(face,
2867  translate_s,
2868  s_sw,
2869  s_ne,
2870  reflected_face,
2871  diff_level,
2872  in_neighbouring_tree);
2873  // Check if they agree...
2874  if (face_neigh_pt != 0)
2875  {
2876  if (face_neigh_pt == edge_neigh_pt)
2877  {
2878  return true;
2879  }
2880  }
2881 
2882  break;
2883 
2884 
2885  case DB:
2886 
2887  // Get first face neighbour
2888  face = D;
2889  face_neigh_pt = gteq_face_neighbour(face,
2890  translate_s,
2891  s_sw,
2892  s_ne,
2893  reflected_face,
2894  diff_level,
2895  in_neighbouring_tree);
2896 
2897  // Check if they agree...
2898  if (face_neigh_pt != 0)
2899  {
2900  if (face_neigh_pt == edge_neigh_pt)
2901  {
2902  return true;
2903  }
2904  }
2905 
2906  // Get second face neighbour
2907  face = B;
2908  face_neigh_pt = gteq_face_neighbour(face,
2909  translate_s,
2910  s_sw,
2911  s_ne,
2912  reflected_face,
2913  diff_level,
2914  in_neighbouring_tree);
2915  // Check if they agree...
2916  if (face_neigh_pt != 0)
2917  {
2918  if (face_neigh_pt == edge_neigh_pt)
2919  {
2920  return true;
2921  }
2922  }
2923 
2924  break;
2925 
2926 
2927  case UB:
2928 
2929  // Get first face neighbour
2930  face = U;
2931  face_neigh_pt = gteq_face_neighbour(face,
2932  translate_s,
2933  s_sw,
2934  s_ne,
2935  reflected_face,
2936  diff_level,
2937  in_neighbouring_tree);
2938 
2939  // Check if they agree...
2940  if (face_neigh_pt != 0)
2941  {
2942  if (face_neigh_pt == edge_neigh_pt)
2943  {
2944  return true;
2945  }
2946  }
2947 
2948  // Get second face neighbour
2949  face = B;
2950  face_neigh_pt = gteq_face_neighbour(face,
2951  translate_s,
2952  s_sw,
2953  s_ne,
2954  reflected_face,
2955  diff_level,
2956  in_neighbouring_tree);
2957 
2958  // Check if they agree...
2959  if (face_neigh_pt != 0)
2960  {
2961  if (face_neigh_pt == edge_neigh_pt)
2962  {
2963  return true;
2964  }
2965  }
2966 
2967  break;
2968 
2969  case LD:
2970 
2971 
2972  // Get first face neighbour
2973  face = L;
2974  face_neigh_pt = gteq_face_neighbour(face,
2975  translate_s,
2976  s_sw,
2977  s_ne,
2978  reflected_face,
2979  diff_level,
2980  in_neighbouring_tree);
2981 
2982  // Check if they agree...
2983  if (face_neigh_pt != 0)
2984  {
2985  if (face_neigh_pt == edge_neigh_pt)
2986  {
2987  return true;
2988  }
2989  }
2990 
2991  // Get second face neighbour
2992  face = D;
2993  face_neigh_pt = gteq_face_neighbour(face,
2994  translate_s,
2995  s_sw,
2996  s_ne,
2997  reflected_face,
2998  diff_level,
2999  in_neighbouring_tree);
3000  // Check if they agree...
3001  if (face_neigh_pt != 0)
3002  {
3003  if (face_neigh_pt == edge_neigh_pt)
3004  {
3005  return true;
3006  }
3007  }
3008 
3009  break;
3010 
3011  case RD:
3012 
3013 
3014  // Get first face neighbour
3015  face = R;
3016  face_neigh_pt = gteq_face_neighbour(face,
3017  translate_s,
3018  s_sw,
3019  s_ne,
3020  reflected_face,
3021  diff_level,
3022  in_neighbouring_tree);
3023 
3024  // Check if they agree...
3025  if (face_neigh_pt != 0)
3026  {
3027  if (face_neigh_pt == edge_neigh_pt)
3028  {
3029  return true;
3030  }
3031  }
3032 
3033  // Get second face neighbour
3034  face = D;
3035  face_neigh_pt = gteq_face_neighbour(face,
3036  translate_s,
3037  s_sw,
3038  s_ne,
3039  reflected_face,
3040  diff_level,
3041  in_neighbouring_tree);
3042  // Check if they agree...
3043  if (face_neigh_pt != 0)
3044  {
3045  if (face_neigh_pt == edge_neigh_pt)
3046  {
3047  return true;
3048  }
3049  }
3050 
3051  break;
3052 
3053  case LU:
3054 
3055  // Get first face neighbour
3056  face = L;
3057  face_neigh_pt = gteq_face_neighbour(face,
3058  translate_s,
3059  s_sw,
3060  s_ne,
3061  reflected_face,
3062  diff_level,
3063  in_neighbouring_tree);
3064 
3065  // Check if they agree...
3066  if (face_neigh_pt != 0)
3067  {
3068  if (face_neigh_pt == edge_neigh_pt)
3069  {
3070  return true;
3071  }
3072  }
3073 
3074  // Get second face neighbour
3075  face = U;
3076  face_neigh_pt = gteq_face_neighbour(face,
3077  translate_s,
3078  s_sw,
3079  s_ne,
3080  reflected_face,
3081  diff_level,
3082  in_neighbouring_tree);
3083 
3084  // Check if they agree...
3085  if (face_neigh_pt != 0)
3086  {
3087  if (face_neigh_pt == edge_neigh_pt)
3088  {
3089  return true;
3090  }
3091  }
3092 
3093  break;
3094 
3095 
3096  case RU:
3097 
3098  // Get first face neighbour
3099  face = R;
3100  face_neigh_pt = gteq_face_neighbour(face,
3101  translate_s,
3102  s_sw,
3103  s_ne,
3104  reflected_face,
3105  diff_level,
3106  in_neighbouring_tree);
3107 
3108  // Check if they agree...
3109  if (face_neigh_pt != 0)
3110  {
3111  if (face_neigh_pt == edge_neigh_pt)
3112  {
3113  return true;
3114  }
3115  }
3116 
3117  // Get second face neighbour
3118  face = U;
3119  face_neigh_pt = gteq_face_neighbour(face,
3120  translate_s,
3121  s_sw,
3122  s_ne,
3123  reflected_face,
3124  diff_level,
3125  in_neighbouring_tree);
3126 
3127  // Check if they agree...
3128  if (face_neigh_pt != 0)
3129  {
3130  if (face_neigh_pt == edge_neigh_pt)
3131  {
3132  return true;
3133  }
3134  }
3135 
3136  break;
3137 
3138 
3139  case LF:
3140 
3141 
3142  // Get first face neighbour
3143  face = L;
3144  face_neigh_pt = gteq_face_neighbour(face,
3145  translate_s,
3146  s_sw,
3147  s_ne,
3148  reflected_face,
3149  diff_level,
3150  in_neighbouring_tree);
3151 
3152  // Check if they agree...
3153  if (face_neigh_pt != 0)
3154  {
3155  if (face_neigh_pt == edge_neigh_pt)
3156  {
3157  return true;
3158  }
3159  }
3160 
3161  // Get second face neighbour
3162  face = F;
3163  face_neigh_pt = gteq_face_neighbour(face,
3164  translate_s,
3165  s_sw,
3166  s_ne,
3167  reflected_face,
3168  diff_level,
3169  in_neighbouring_tree);
3170 
3171  // Check if they agree...
3172  if (face_neigh_pt != 0)
3173  {
3174  if (face_neigh_pt == edge_neigh_pt)
3175  {
3176  return true;
3177  }
3178  }
3179 
3180  break;
3181 
3182  case RF:
3183 
3184  // Get first face neighbour
3185  face = R;
3186  face_neigh_pt = gteq_face_neighbour(face,
3187  translate_s,
3188  s_sw,
3189  s_ne,
3190  reflected_face,
3191  diff_level,
3192  in_neighbouring_tree);
3193 
3194  // Check if they agree...
3195  if (face_neigh_pt != 0)
3196  {
3197  if (face_neigh_pt == edge_neigh_pt)
3198  {
3199  return true;
3200  }
3201  }
3202 
3203  // Get second face neighbour
3204  face = F;
3205  face_neigh_pt = gteq_face_neighbour(face,
3206  translate_s,
3207  s_sw,
3208  s_ne,
3209  reflected_face,
3210  diff_level,
3211  in_neighbouring_tree);
3212 
3213  // Check if they agree...
3214  if (face_neigh_pt != 0)
3215  {
3216  if (face_neigh_pt == edge_neigh_pt)
3217  {
3218  return true;
3219  }
3220  }
3221 
3222  break;
3223 
3224 
3225  case DF:
3226 
3227  // Get first face neighbour
3228  face = D;
3229  face_neigh_pt = gteq_face_neighbour(face,
3230  translate_s,
3231  s_sw,
3232  s_ne,
3233  reflected_face,
3234  diff_level,
3235  in_neighbouring_tree);
3236 
3237  // Check if they agree...
3238  if (face_neigh_pt != 0)
3239  {
3240  if (face_neigh_pt == edge_neigh_pt)
3241  {
3242  return true;
3243  }
3244  }
3245 
3246 
3247  // Get second face neighbour
3248  face = F;
3249  face_neigh_pt = gteq_face_neighbour(face,
3250  translate_s,
3251  s_sw,
3252  s_ne,
3253  reflected_face,
3254  diff_level,
3255  in_neighbouring_tree);
3256 
3257  // Check if they agree...
3258  if (face_neigh_pt != 0)
3259  {
3260  if (face_neigh_pt == edge_neigh_pt)
3261  {
3262  return true;
3263  }
3264  }
3265 
3266  break;
3267 
3268 
3269  case UF:
3270 
3271  // Get first face neighbour
3272  face = U;
3273  face_neigh_pt = gteq_face_neighbour(face,
3274  translate_s,
3275  s_sw,
3276  s_ne,
3277  reflected_face,
3278  diff_level,
3279  in_neighbouring_tree);
3280 
3281  // Check if they agree...
3282  if (face_neigh_pt != 0)
3283  {
3284  if (face_neigh_pt == edge_neigh_pt)
3285  {
3286  return true;
3287  }
3288  }
3289 
3290 
3291  // Get second face neighbour
3292  face = F;
3293  face_neigh_pt = gteq_face_neighbour(face,
3294  translate_s,
3295  s_sw,
3296  s_ne,
3297  reflected_face,
3298  diff_level,
3299  in_neighbouring_tree);
3300 
3301  // Check if they agree...
3302  if (face_neigh_pt != 0)
3303  {
3304  if (face_neigh_pt == edge_neigh_pt)
3305  {
3306  return true;
3307  }
3308  }
3309 
3310  break;
3311 
3312  default:
3313 
3314  // There is no face neighbour in this direction so they can't
3315  // agree:
3316  std::ostringstream error_stream;
3317  error_stream << "Never get here! Edge:" << Direct_string[edge] << " "
3318  << edge << std::endl;
3319  throw OomphLibError(
3320  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3321  }
3322 
3323  // If we've made it to here then we've located the requested edge
3324  // but found that none of its two face neighbours match the specified
3325  // edge neighbour:
3326  return false;
3327  }
D
dominoes D
Definition: Domino.cpp:55
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:48
oomph::OcTree::gteq_face_neighbour
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
Definition: octree.cc:3373
oomph::OcTreeNames::DF
@ DF
Definition: octree.h:67
oomph::OcTreeNames::RD
@ RD
Definition: octree.h:62
oomph::OcTreeNames::RB
@ RB
Definition: octree.h:58
oomph::OcTreeNames::RU
@ RU
Definition: octree.h:64
oomph::OcTreeNames::LF
@ LF
Definition: octree.h:65
oomph::OcTreeNames::LD
@ LD
Definition: octree.h:61
oomph::OcTreeNames::RF
@ RF
Definition: octree.h:66
oomph::OcTreeNames::LU
@ LU
Definition: octree.h:63
oomph::OcTreeNames::DB
@ DB
Definition: octree.h:59
oomph::OcTreeNames::UB
@ UB
Definition: octree.h:60

References D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, Direct_string, oomph::OcTreeNames::F, gteq_face_neighbour(), L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, and oomph::OcTreeNames::UF.

Referenced by gteq_true_edge_neighbour().

◆ faces_of_common_edge()

Vector< int > oomph::OcTree::faces_of_common_edge ( const int &  edge)
static

Function that, given an edge, returns the two faces on which it.

Given an edge, this function returns the faces on which it lies.

269  {
270  using namespace OcTreeNames;
271 
272 #ifdef PARANOID
273  if ((edge != LB) && (edge != RB) && (edge != DB) && (edge != UB) &&
274  (edge != LD) && (edge != RD) && (edge != LU) && (edge != RU) &&
275  (edge != LF) && (edge != RF) && (edge != DF) && (edge != UF))
276  {
277  std::ostringstream error_stream;
278  error_stream << "Edge " << Direct_string[edge] << "is not valid!";
279  throw OomphLibError(
280  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
281  }
282 #endif
283 
284  // Allocate space for the result
285  Vector<int> faces(2, 0);
286 
287  // Which faces does this edge lie on?
288  switch (edge)
289  {
290  case LB:
291  // First entry
292  faces[0] = L;
293 
294  // Second entry
295  faces[1] = B;
296 
297  // Break
298  break;
299  case RB:
300  // First entry
301  faces[0] = R;
302 
303  // Second entry
304  faces[1] = B;
305 
306  // Break
307  break;
308  case DB:
309  // First entry
310  faces[0] = D;
311 
312  // Second entry
313  faces[1] = B;
314 
315  // Break
316  break;
317  case UB:
318  // First entry
319  faces[0] = U;
320 
321  // Second entry
322  faces[1] = B;
323 
324  // Break
325  break;
326  case LD:
327  // First entry
328  faces[0] = L;
329 
330  // Second entry
331  faces[1] = D;
332 
333  // Break
334  break;
335  case RD:
336  // First entry
337  faces[0] = R;
338 
339  // Second entry
340  faces[1] = D;
341 
342  // Break
343  break;
344  case LU:
345  // First entry
346  faces[0] = L;
347 
348  // Second entry
349  faces[1] = U;
350 
351  // Break
352  break;
353  case RU:
354  // First entry
355  faces[0] = R;
356 
357  // Second entry
358  faces[1] = U;
359 
360  // Break
361  break;
362  case LF:
363  // First entry
364  faces[0] = L;
365 
366  // Second entry
367  faces[1] = F;
368 
369  // Break
370  break;
371  case RF:
372  // First entry
373  faces[0] = R;
374 
375  // Second entry
376  faces[1] = F;
377 
378  // Break
379  break;
380  case DF:
381  // First entry
382  faces[0] = D;
383 
384  // Second entry
385  faces[1] = F;
386 
387  // Break
388  break;
389  case UF:
390  // First entry
391  faces[0] = U;
392 
393  // Second entry
394  faces[1] = F;
395 
396  // Break
397  break;
398  default:
399  // Throw an error
400  throw OomphLibError("Incorrect edge input!",
401  OOMPH_CURRENT_FUNCTION,
402  OOMPH_EXCEPTION_LOCATION);
403  }
404 
405  // Return the faces
406  return faces;
407  } // End of faces_of_common_edge

References D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, Direct_string, oomph::OcTreeNames::F, L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, and oomph::OcTreeNames::UF.

Referenced by doc_true_edge_neighbours(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::node_created_by_neighbour(), oomph::RefineableQElement< 3 >::node_created_by_neighbour(), and oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::node_created_by_son_of_neighbour().

◆ get_the_other_face()

int oomph::OcTree::get_the_other_face ( const unsigned &  n1,
const unsigned &  n2,
const unsigned &  nnode1d,
const int &  face 
)
static

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 face and returns the other one. (face is a direction in the set U,D,F,B,L,R).

This function takes as argument two node numbers of two nodes delimiting an edge, and one face of this edge and returns the other face that is sharing this edge. The node numbers given to this function MUST be vertices nodes to work. it also need the value of nnode1d to work. (face is a direction in the set U,D,F,B,L,R).

569  {
570  Vector<int> vect_node1(3);
571  Vector<int> vect_node2(3);
572  Vector<int> vect_face(3);
573  Vector<int> vect_other_face(3);
574 
575  // Translate the nodes to vectors
576  vect_node1 = vertex_node_to_vector(n1, nnode1d);
577  vect_node2 = vertex_node_to_vector(n2, nnode1d);
578 
579  // Translate the face to a vector
580  vect_face = Direction_to_vector[face];
581 
582  // Compute the vector to the other face -- magic, courtesy of Renaud
583  // Schleck!
584  for (unsigned i = 0; i < 3; i++)
585  {
586  // Calculate the i-th entry
587  vect_other_face[i] = (vect_node1[i] + vect_node2[i]) / 2 - vect_face[i];
588 
589 #ifdef PARANOID
590  if ((vect_other_face[i] != 1) && (vect_other_face[i] != -1) &&
591  (vect_other_face[i] != 0))
592  {
593  throw OomphLibError("The nodes given are not vertices",
594  OOMPH_CURRENT_FUNCTION,
595  OOMPH_EXCEPTION_LOCATION);
596  }
597 #endif
598  }
599 
600  // return the corresponding face
601  return Vector_to_direction[vect_other_face];
602  }
oomph::OcTree::vertex_node_to_vector
static Vector< int > vertex_node_to_vector(const unsigned &n, const unsigned &nnode1d)
Definition: octree.cc:232
oomph::OcTree::Direction_to_vector
static Vector< Vector< int > > Direction_to_vector
Definition: octree.h:353
oomph::OcTree::Vector_to_direction
static std::map< Vector< int >, int > Vector_to_direction
Definition: octree.h:348

References Direction_to_vector, i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Vector_to_direction, and vertex_node_to_vector().

Referenced by oomph::OcTreeForest::construct_up_right_equivalents().

◆ gteq_edge_neighbour()

OcTree * oomph::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
private

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). i_root_edge_neighbour specifies which of these is used. Use this as "reverse communication": First call with i_root_edge_neighbour=0 and n_root_edge_neighour initialised to anything you want (zero, ideally). On return from the fct, n_root_edge_neighour contains the total number of true edge neighbours, so additional calls to the fct with 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.

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!

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 previous son types and the maximum level to which search is performed.

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). i_root_edge_neighbour specifies which of these is used. Use this as "reverse communication": First call with i_root_edge_neighbour=0 and n_root_edge_neighour initialised to anything you want (zero, ideally). On return from the fct, n_root_edge_neighour contains the total number of true edge neighbours, so additional calls to the fct with i_root_edge_neighbour>0 can be made until they've all been visited.
  • s_diff: 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.
4022  {
4023  using namespace OcTreeNames;
4024 
4025 
4026 #ifdef PARANOID
4027  if ((direction != LB) && (direction != RB) && (direction != DB) &&
4028  (direction != UB) && (direction != LD) && (direction != RD) &&
4029  (direction != LU) && (direction != RU) && (direction != LF) &&
4030  (direction != RF) && (direction != DF) && (direction != UF))
4031  {
4032  std::ostringstream error_stream;
4033  error_stream << "Wrong direction: " << Direct_string[direction]
4034  << std::endl;
4035  throw OomphLibError(
4036  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4037  }
4038 #endif
4039 
4040  // Initialise total number of edge neighbours available across
4041  // edges of octree roots
4042  nroot_edge_neighbour = 0;
4043 
4044  OcTree* next_el_pt = 0;
4045  OcTree* return_el_pt = 0;
4046 
4047  // STEP 1: Find the common ancestor
4048  //--------
4049  // Does the element have a father?
4050  if (Father_pt != 0)
4051  {
4052  // If the present octant (whose location inside its
4053  // father element is specified by Son_type) is adjacent to the
4054  // father's edge in the required direction, then its neighbour has
4055  // a different father ---> we need to climb up the tree to
4056  // the father and find his neighbour in the required direction
4057  if (Is_adjacent(direction, Son_type))
4058  {
4059  next_el_pt = dynamic_cast<OcTree*>(Father_pt)->gteq_edge_neighbour(
4060  direction,
4061  i_root_edge_neighbour,
4062  nroot_edge_neighbour,
4063  s_diff,
4064  diff_level,
4065  max_level,
4066  orig_root_pt);
4067  }
4068  // If the present octant (whose location inside its
4069  // father element is specified by Son_type) is adjacent to the
4070  // father's face in the required direction, then its neighbour has
4071  // a different father ---> we need to climb up the tree to
4072  // the father and find his neighbour in the required direction,
4073  // crossing the face as we do so.
4074  else if (Common_face(direction, Son_type) != OMEGA)
4075  {
4076  // Initialise bool
4077  bool in_neighbouring_tree = false;
4078 
4079  // We're going across a face:
4080  double s_difflo = 0.0;
4081  double s_diffhi = 0.0;
4082  int diff_level_edge = 0;
4083 
4084  next_el_pt = dynamic_cast<OcTree*>(Father_pt)->gteq_face_neighbour(
4085  Common_face(direction, Son_type),
4086  s_difflo,
4087  s_diffhi,
4088  diff_level_edge,
4089  in_neighbouring_tree,
4090  max_level,
4091  orig_root_pt);
4092  }
4093  // If the present octant is not adjacent to the
4094  // father's face/edge in the required direction, then the
4095  // neighbour has the same father and will be obtained
4096  // by the appropriate reflection inside the father element
4097  else
4098  {
4099  next_el_pt = dynamic_cast<OcTree*>(Father_pt);
4100  }
4101 
4102  // We're about to ascend one level:
4103  diff_level -= 1;
4104 
4105  // Work out position of "low" vertex of present edge
4106  // in its father element
4107  s_diff += pow(0.5, -diff_level) * S_direct_edge(direction, Son_type);
4108 
4109  // STEP 2: We have now located the neighbour's father and need to
4110  // ------- find the appropriate son.
4111  // Buffer cases where ...
4112  if (next_el_pt != 0)
4113  {
4114  // If the father is a leaf then we can't descend to the same
4115  // level as the present node ---> simply return the father himself
4116  // as the (greater) neighbour. Same applies if we are about
4117  // to descend lower than the max_level (in a neighbouring tree)
4118  if ((next_el_pt->Son_pt.size() == 0) ||
4119  (next_el_pt->Level > max_level - 1))
4120  {
4121  return_el_pt = next_el_pt;
4122  }
4123  // We have located the neighbour's father: The position of the
4124  // neighbour is obtained by `reflecting' the position of the
4125  // node itself.
4126  else
4127  {
4128  // By default (in the absence of rotations) we obtain the
4129  // son octant by reflecting
4130  int son_octant = Reflect(direction, Son_type);
4131 
4132  // If there is a rotation, we rotate the son octant
4133  if (orig_root_pt != next_el_pt->Root_pt)
4134  {
4135  int my_up = dynamic_cast<OcTreeRoot*>(Root_pt)->up_equivalent(
4136  next_el_pt->Root_pt);
4137  int my_right = dynamic_cast<OcTreeRoot*>(Root_pt)->right_equivalent(
4138  next_el_pt->Root_pt);
4139  son_octant = rotate(my_up, my_right, son_octant);
4140  }
4141 
4142  return_el_pt = dynamic_cast<OcTree*>(next_el_pt->Son_pt[son_octant]);
4143 
4144  // Work out position of "low" vertex of present edge
4145  // in next higher element
4146  s_diff -= pow(0.5, -diff_level) * S_direct_edge(direction, Son_type);
4147 
4148  // We have just descended one level
4149  diff_level += 1;
4150  }
4151  }
4152  // The neighbour's father lies outside the boundary --> the neighbour
4153  // itself does too --> return NULL.
4154  else
4155  {
4156  return_el_pt = 0;
4157  }
4158  }
4159  // Element does not have a father --> check if it has a neighbouring
4160  // tree in the appropriate direction
4161  else
4162  {
4163  // Get total number of edge neighbours available across
4164  // edges of octree roots for return
4165  nroot_edge_neighbour =
4166  dynamic_cast<OcTreeRoot*>(Root_pt)->nedge_neighbour(direction);
4167 
4168  // Get vector of edge neighbours (if any) in appropriate direction
4169  Vector<TreeRoot*> edge_neighbour_pt =
4170  dynamic_cast<OcTreeRoot*>(Root_pt)->edge_neighbour_pt(direction);
4171 
4172  // Get the number of edge neighbours
4173  unsigned n_neigh = edge_neighbour_pt.size();
4174 
4175  // If there are any edge neighbours
4176  if (n_neigh > 0)
4177  {
4178  // Return the appropriate edge neighbour
4179  return_el_pt =
4180  dynamic_cast<OcTree*>(edge_neighbour_pt[i_root_edge_neighbour]);
4181  }
4182  else
4183  {
4184  return_el_pt = 0;
4185  }
4186  }
4187 
4188  return return_el_pt;
4189  } // End of gteq_edge_neighbour
oomph::OcTree::Is_adjacent
static DenseMatrix< bool > Is_adjacent
Definition: octree.h:529
oomph::OcTree::rotate
static Vector< int > rotate(const int &new_up, const int &new_right, const Vector< int > &dir)
Definition: octree.cc:750
oomph::OcTree::Common_face
static DenseMatrix< int > Common_face
Definition: octree.h:537
oomph::OcTree::Reflect
static DenseMatrix< int > Reflect
Definition: octree.h:533
oomph::OcTree::gteq_edge_neighbour
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
Definition: octree.cc:4015
oomph::OcTree::S_direct_edge
static DenseMatrix< double > S_direct_edge
Definition: octree.h:600
oomph::Tree::Father_pt
Tree * Father_pt
Pointer to the Father of the Tree.
Definition: tree.h:296
oomph::Tree::Root_pt
TreeRoot * Root_pt
Pointer to the root of the tree.
Definition: tree.h:292
oomph::Tree::Son_type
int Son_type
Son type (e.g. SW/SE/NW/NE in a quadtree)
Definition: tree.h:305

References Common_face, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, Direct_string, oomph::Tree::Father_pt, gteq_face_neighbour(), Is_adjacent, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::Tree::Level, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, oomph::Tree::OMEGA, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, Reflect, oomph::OcTreeNames::RF, oomph::Tree::Root_pt, rotate(), oomph::OcTreeNames::RU, S_direct_edge, oomph::Tree::Son_pt, oomph::Tree::Son_type, oomph::OcTreeNames::UB, and oomph::OcTreeNames::UF.

Referenced by gteq_true_edge_neighbour().

◆ gteq_face_neighbour() [1/2]

OcTree * oomph::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
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.

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 previous son types and the maximum level to which search is performed.

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.
  • face: We're looking for the neighbour across our face 'direction' (L/R/U/D/B/F). When viewed from the neighbour, this face is ‘face’ (L/R/U/D/B/F). [If there's no relative rotation between neighbours then this is a mere reflection, e.g. direction=F --> face=B etc.]
  • 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.
3846  {
3847  using namespace OcTreeNames;
3848 
3849 #ifdef PARANOID
3850  if ((direction != L) && (direction != R) && (direction != F) &&
3851  (direction != B) && (direction != U) && (direction != D))
3852  {
3853  std::ostringstream error_stream;
3854  error_stream << "Direction " << Direct_string[direction]
3855  << " is not L, R, B, F, D or U." << std::endl;
3856  throw OomphLibError(
3857  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3858  }
3859 #endif
3860 
3861  OcTree* next_el_pt;
3862  OcTree* return_el_pt;
3863 
3864  // Initialise in_neighbouring tree to false. It will be set to true
3865  // during the recursion if we do actually hop over in to the neighbour
3866  in_neighbouring_tree = false;
3867 
3868  // STEP 1: Find the neighbour's father
3869  //--------
3870  // Does the element have a father?
3871  if (Father_pt != 0)
3872  {
3873  // If the present octant (whose location inside its
3874  // father element is specified by Son_type) is adjacent to the
3875  // father's face in the required direction, then its neighbour has
3876  // a different father ---> we need to climb up the tree to
3877  // the father and find his neighbour in the required direction
3878  if (Is_adjacent(direction, Son_type))
3879  {
3880  next_el_pt = dynamic_cast<OcTree*>(Father_pt)->gteq_face_neighbour(
3881  direction,
3882  s_difflo,
3883  s_diffhi,
3884  diff_level,
3885  in_neighbouring_tree,
3886  max_level,
3887  orig_root_pt);
3888  }
3889  // If the present octant is not adjacent to the
3890  // father's face in the required direction, then the
3891  // neighbour has the same father and will be obtained
3892  // by the appropriate reflection inside the father element
3893  else
3894  {
3895  next_el_pt = dynamic_cast<OcTree*>(Father_pt);
3896  }
3897 
3898  // We're about to ascend one level:
3899  diff_level -= 1;
3900 
3901  // Work out position of lower-left corner of present face
3902  // in its father element
3903  s_difflo += pow(0.5, -diff_level) * S_directlo(direction, Son_type);
3904  s_diffhi += pow(0.5, -diff_level) * S_directhi(direction, Son_type);
3905 
3906  // STEP 2: We have now located the neighbour's father and need to
3907  // ------- find the appropriate son.
3908  // Buffer cases where the neighbour (and hence its father) lie outside
3909  // the boundary
3910  if (next_el_pt != 0)
3911  {
3912  // If the father is a leaf then we can't descend to the same
3913  // level as the present node ---> simply return the father himself
3914  // as the (greater) neighbour. Same applies if we are about
3915  // to descend lower than the max_level (in a neighbouring tree)
3916  if ((next_el_pt->Son_pt.size() == 0) ||
3917  (next_el_pt->Level > max_level - 1))
3918  {
3919  return_el_pt = next_el_pt;
3920  }
3921  // We have located the neighbour's father: The position of the
3922  // neighbour is obtained by `reflecting' the position of the
3923  // node itself.
3924  else
3925  {
3926  // By default (in the absence of rotations) we obtain the
3927  // son octant by reflecting
3928  int son_octant = Reflect(direction, Son_type);
3929 
3930  // If there is a rotation, we rotate the son octant
3931  if (orig_root_pt != next_el_pt->Root_pt)
3932  {
3933  int my_up = dynamic_cast<OcTreeRoot*>(Root_pt)->up_equivalent(
3934  next_el_pt->Root_pt);
3935  int my_right = dynamic_cast<OcTreeRoot*>(Root_pt)->right_equivalent(
3936  next_el_pt->Root_pt);
3937  son_octant = rotate(my_up, my_right, son_octant);
3938  }
3939 
3940  return_el_pt = dynamic_cast<OcTree*>(next_el_pt->Son_pt[son_octant]);
3941 
3942  // Work out position of lower-left corner of present face
3943  // in next higher element
3944  s_difflo -= pow(0.5, -diff_level) * S_directlo(direction, Son_type);
3945  s_diffhi -= pow(0.5, -diff_level) * S_directhi(direction, Son_type);
3946 
3947  // We have just descended one level
3948  diff_level += 1;
3949  }
3950  }
3951  // The neighbour's father lies outside the boundary --> the neighbour
3952  // itself does too --> return NULL.
3953  else
3954  {
3955  return_el_pt = 0;
3956  }
3957  }
3958  // Element does not have a father --> check if it has a neighbouring
3959  // tree in the appropriate direction
3960  else
3961  {
3962  // Find neighbouring root
3963  if (Root_pt->neighbour_pt(direction) != 0)
3964  {
3965  // If we're in the neighbouring tree
3966  in_neighbouring_tree = true;
3967 
3968  // Return
3969  return_el_pt = dynamic_cast<OcTree*>(Root_pt->neighbour_pt(direction));
3970  }
3971  // No neighbouring tree, so there really is no neighbour --> return NULL
3972  else
3973  {
3974  return_el_pt = 0;
3975  }
3976  }
3977 
3978  // Return the appropriate OcTree pointer
3979  return return_el_pt;
3980  } // End of gteq_face_neighbour
oomph::Matrix
Definition: matrices.h:74
oomph::OcTree::S_directhi
static DenseMatrix< double > S_directhi
Definition: octree.h:581
oomph::OcTree::S_directlo
static DenseMatrix< double > S_directlo
Definition: octree.h:573
oomph::TreeRoot::neighbour_pt
TreeRoot *& neighbour_pt(const int &direction)
Definition: tree.h:357

References D, Direct_string, oomph::OcTreeNames::F, oomph::Tree::Father_pt, gteq_face_neighbour(), Is_adjacent, L, oomph::Tree::Level, oomph::TreeRoot::neighbour_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), R, Reflect, oomph::Tree::Root_pt, rotate(), S_directhi, S_directlo, oomph::Tree::Son_pt, oomph::Tree::Son_type, and RachelsAdvectionDiffusion::U.

◆ gteq_face_neighbour() [2/2]

OcTree * oomph::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 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 General convention:

  • 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 translate_s turns the index of the local coordinate in the present octree into that of the neighbour. If there are no rotations then 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 (s_sw[0], s_sw[1], 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 (s_ne[0], s_ne[1], s_ne[2]) in the neighbouring octree.
  • We're looking for a neighbour in the specified direction. When viewed from the neighbouring octree, the face that separates the present octree from its neighbour is the neighbour's face 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 R [ight] direction, face will be L [eft]
  • diff_level <= 0 indicates the difference in refinement levels between the two neighbours. If diff_level==0, the neighbour has the same size as the current octree.
3380  {
3381  using namespace OcTreeNames;
3382 
3383 #ifdef PARANOID
3384  if ((direction != L) && (direction != R) && (direction != F) &&
3385  (direction != B) && (direction != U) && (direction != D))
3386  {
3387  std::ostringstream error_stream;
3388  error_stream << "Wrong direction: " << Direct_string[direction]
3389  << std::endl;
3390  throw OomphLibError(
3391  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3392  }
3393 #endif
3394 
3395  // Initialise in_neighbouring tree to false. It will be set to true
3396  // during the recursion if we do actually hop over in to the neighbour
3397  in_neighbouring_tree = false;
3398 
3399  // Maximum level to which we're allowed to descend (we only want
3400  // greater-or-equal-sized neighbours)
3401  int max_level = Level;
3402 
3403  // Current element has the following root:
3404  OcTreeRoot* orig_root_pt = dynamic_cast<OcTreeRoot*>(Root_pt);
3405 
3406  // Initialise offset in local coordinate
3407  double s_difflo = 0;
3408  double s_diffhi = 0;
3409 
3410  // Initialise difference in level
3411  diff_level = 0;
3412 
3413  // Find neighbour
3414  OcTree* return_pt = gteq_face_neighbour(direction,
3415  s_difflo,
3416  s_diffhi,
3417  diff_level,
3418  in_neighbouring_tree,
3419  max_level,
3420  orig_root_pt);
3421 
3422  OcTree* neighb_pt = return_pt;
3423 
3424  // Initialise the translation scheme
3425  for (unsigned i = 0; i < 3; i++)
3426  {
3427  translate_s[i] = i;
3428  }
3429 
3430  // If neighbour exists: What's the direction of the interfacial
3431  // face when viewed from within the neighbour element?
3432  if (neighb_pt != 0)
3433  {
3434  // Find the reflection of the original direction, which will be the
3435  // direction to the face in the neighbour, if there are no rotations.
3436  int reflected_dir = Reflect_face[direction];
3437 
3438  // These coordinates are the coordinates of the ne and sw points
3439  // in the neighbour (provided there are no rotations -- we'll correct
3440  // for these below)
3441  s_sw[0] = S_base(0, reflected_dir) +
3442  S_steplo(0, reflected_dir) * s_difflo +
3443  S_stephi(0, reflected_dir) * s_diffhi;
3444 
3445  s_sw[1] = S_base(1, reflected_dir) +
3446  S_steplo(1, reflected_dir) * s_difflo +
3447  S_stephi(1, reflected_dir) * s_diffhi;
3448 
3449  s_sw[2] = S_base(2, reflected_dir) +
3450  S_steplo(2, reflected_dir) * s_difflo +
3451  S_stephi(2, reflected_dir) * s_diffhi;
3452 
3453  s_ne[0] = S_base(0, reflected_dir) +
3454  S_steplo(0, reflected_dir) * pow(2.0, diff_level) +
3455  S_steplo(0, reflected_dir) * s_difflo +
3456  S_stephi(0, reflected_dir) * pow(2.0, diff_level) +
3457  S_stephi(0, reflected_dir) * s_diffhi;
3458 
3459  s_ne[1] = S_base(1, reflected_dir) +
3460  S_steplo(1, reflected_dir) * pow(2.0, diff_level) +
3461  S_steplo(1, reflected_dir) * s_difflo +
3462  S_stephi(1, reflected_dir) * pow(2.0, diff_level) +
3463  S_stephi(1, reflected_dir) * s_diffhi;
3464 
3465  s_ne[2] = S_base(2, reflected_dir) +
3466  S_steplo(2, reflected_dir) * pow(2.0, diff_level) +
3467  S_steplo(2, reflected_dir) * s_difflo +
3468  S_stephi(2, reflected_dir) * pow(2.0, diff_level) +
3469  S_stephi(2, reflected_dir) * s_diffhi;
3470 
3471  // If there is no rotation then the new direction is the same as the
3472  // old direction
3473  int new_dir = direction;
3474 
3475  // If necessary, rotate things around (the orientation of Up in the
3476  // neighbour might be be different from that in the present element)
3477  // If the root of the neighbour is the not same as the one of the present
3478  // element then their orientations may not be the same and the new
3479  // direction is given by :
3480  if (neighb_pt->Root_pt != Root_pt)
3481  {
3482  new_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3483  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3484  direction);
3485  }
3486 
3487  // What's the direction of the interfacial face when viewed from within
3488  // the neighbour element?
3489  face = Reflect_face[new_dir];
3490 
3491  Vector<double> s_sw_new(3), s_ne_new(3);
3492 
3493  // If the root of the present element is different from the root
3494  // of his neighbour, we have to rotate the RUF and LDB coordinates
3495  // to have their equivalents in the neighbour's point of view.
3496  if (neighb_pt->Root_pt != Root_pt)
3497  {
3498  int tmp_dir;
3499  Vector<int> vect1(3);
3500  Vector<int> vect2(3);
3501  Vector<int> vect3(3);
3502  DenseMatrix<int> Mat_rot(3);
3503 
3504  // All this is just to compute the rotation matrix
3505  tmp_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3506  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3507  R);
3508  vect1 = Direction_to_vector[tmp_dir];
3509 
3510  // All this is just to compute the rotation matrix
3511  tmp_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3512  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3513  U);
3514  vect2 = Direction_to_vector[tmp_dir];
3515 
3516  // All this is just to compute the rotation matrix
3517  tmp_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3518  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3519  F);
3520  vect3 = Direction_to_vector[tmp_dir];
3521 
3522  // Setup the inverse rotation matrix
3523  for (int i = 0; i < 3; i++)
3524  {
3525  Mat_rot(i, 0) = vect1[i];
3526  Mat_rot(i, 1) = vect2[i];
3527  Mat_rot(i, 2) = vect3[i];
3528  }
3529 
3530  // Initialise the translation scheme
3531  Vector<int> translate_s_new(3);
3532 
3533  // Then the rotation of the coordinates
3534  for (unsigned i = 0; i < 3; i++)
3535  {
3536  s_ne_new[i] = 0.0;
3537  s_sw_new[i] = 0.0;
3538  translate_s_new[i] = 0;
3539  for (unsigned k = 0; k < 3; k++)
3540  {
3541  s_ne_new[i] += Mat_rot(i, k) * s_ne[k];
3542  s_sw_new[i] += Mat_rot(i, k) * s_sw[k];
3543  translate_s_new[i] += Mat_rot(i, k) * translate_s[k];
3544  }
3545  }
3546 
3547  s_ne = s_ne_new;
3548  s_sw = s_sw_new;
3549 
3550  // Set the translation scheme
3551  for (unsigned i = 0; i < 3; i++)
3552  {
3553  // abs is ok here; not fabs!
3554  translate_s[i] = std::abs(translate_s_new[i]);
3555  }
3556  }
3557 
3558  } // end of if(neighb_pt!=0)
3559 
3560  return return_pt;
3561  }
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
oomph::OcTree::Reflect_face
static Vector< int > Reflect_face
Get opposite face, e.g. Reflect_face[L]=R.
Definition: octree.h:332
oomph::Tree::Level
int Level
Level of the Tree (level 0 = root)
Definition: tree.h:302
k
char char char int int * k
Definition: level2_impl.h:374

References abs(), D, Direct_string, Direction_to_vector, oomph::OcTreeNames::F, i, k, L, oomph::Tree::Level, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), R, Reflect_face, oomph::OcTreeRoot::right_equivalent(), oomph::Tree::Root_pt, rotate(), S_base, S_stephi, S_steplo, RachelsAdvectionDiffusion::U, and oomph::OcTreeRoot::up_equivalent().

Referenced by oomph::RefineableQElement< 3 >::check_integrity(), doc_face_neighbours(), edge_neighbour_is_face_neighbour(), gteq_edge_neighbour(), gteq_face_neighbour(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::node_created_by_neighbour(), oomph::RefineableQElement< 3 >::node_created_by_neighbour(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::node_created_by_son_of_neighbour(), oomph::RefineableQElement< 3 >::oc_hang_helper(), and oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::oc_hang_helper().

◆ gteq_true_edge_neighbour()

OcTree * oomph::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

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). i_root_edge_neighbour specifies which of these is used. Use this as "reverse communication": First call with i_root_edge_neighbour=0 and n_root_edge_neighour initialised to anything you want (zero, ideally). On return from the fct, n_root_edge_neighour contains the total number of true edge neighbours, so additional calls to the fct with i_root_edge_neighbour>0 can be made until they've all been visited.
  • The vector translate_s turns the index of the local coordinate in the present octree into that of the neighbour. If there are no rotations then 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 (s_lo[0], s_lo[1], 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 (s_hi[0], s_hi[1], s_hi[2])
  • We're looking for a neighbour in the specified direction. When viewed from the neighbouring octree, the edge that separates the present octree from its neighbour is the neighbour's edge 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 DB direction, edge will be UF.
  • diff_level <= 0 indicates the difference in refinement levels between the two neighbours. If diff_level==0, the neighbour has the same size as the current octree.

Important: We're only looking for 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.

3627  {
3628  using namespace OcTreeNames;
3629 
3630 #ifdef PARANOID
3631  if ((direction != LB) && (direction != RB) && (direction != DB) &&
3632  (direction != UB) && (direction != LD) && (direction != RD) &&
3633  (direction != LU) && (direction != RU) && (direction != LF) &&
3634  (direction != RF) && (direction != DF) && (direction != UF))
3635  {
3636  std::ostringstream error_stream;
3637  error_stream << "Wrong direction: " << Direct_string[direction]
3638  << std::endl;
3639  throw OomphLibError(
3640  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3641  }
3642 #endif
3643 
3644  // Maximum level to which we're allowed to descend (we only want
3645  // greater-or-equal-sized neighbours)
3646  int max_level = Level;
3647 
3648  // Current element has the following root:
3649  OcTreeRoot* orig_root_pt = dynamic_cast<OcTreeRoot*>(Root_pt);
3650 
3651  // Initialise offset in local coordinate along edge
3652  double s_diff = 0;
3653 
3654  // Initialise difference in level
3655  diff_level = 0;
3656 
3657  // Find edge neighbour
3658  OcTree* return_pt = gteq_edge_neighbour(direction,
3659  i_root_edge_neighbour,
3660  nroot_edge_neighbour,
3661  s_diff,
3662  diff_level,
3663  max_level,
3664  orig_root_pt);
3665 
3666  // Only use "true" edge neighbours
3667  if (edge_neighbour_is_face_neighbour(direction, return_pt))
3668  {
3669  return_pt = 0;
3670  }
3671 
3672  // By default, we return what was returned as the true edge neighbour.
3673  OcTree* neighb_pt = return_pt;
3674 
3675  // Initialise the translation scheme
3676  for (unsigned i = 0; i < 3; i++)
3677  {
3678  translate_s[i] = i;
3679  }
3680 
3681  // If neighbour exists: What's the direction of the interfacial
3682  // edge when viewed from within the neighbour element?
3683  if (neighb_pt != 0)
3684  {
3685  // Find the reflection of the original direction, which will be the
3686  // direction to the edge in the neighbour, if there are no rotations.
3687  int reflected_dir = Reflect_edge[direction];
3688 
3689  // These coordinates are the coordinates of the "low" and "high" points
3690  // in the neighbour (provided there are no rotations -- we'll correct
3691  // for these below)
3692  s_lo[0] =
3693  S_base_edge(0, reflected_dir) + S_step_edge(0, reflected_dir) * s_diff;
3694 
3695  s_lo[1] =
3696  S_base_edge(1, reflected_dir) + S_step_edge(1, reflected_dir) * s_diff;
3697 
3698  s_lo[2] =
3699  S_base_edge(2, reflected_dir) + S_step_edge(2, reflected_dir) * s_diff;
3700 
3701  s_hi[0] = S_base_edge(0, reflected_dir) +
3702  S_step_edge(0, reflected_dir) * pow(2.0, diff_level) +
3703  S_step_edge(0, reflected_dir) * s_diff;
3704 
3705  s_hi[1] = S_base_edge(1, reflected_dir) +
3706  S_step_edge(1, reflected_dir) * pow(2.0, diff_level) +
3707  S_step_edge(1, reflected_dir) * s_diff;
3708 
3709  s_hi[2] = S_base_edge(2, reflected_dir) +
3710  S_step_edge(2, reflected_dir) * pow(2.0, diff_level) +
3711  S_step_edge(2, reflected_dir) * s_diff;
3712 
3713  // If there is no rotation then the new direction is the same as the
3714  // old direction
3715  int new_dir = direction;
3716 
3717 
3718  // If necessary, rotate things around (the orientation of Up in the
3719  // neighbour might be be different from that in the present element)
3720  // If the root of the neighbour is the not same as the one of the present
3721  // element then their orientations may not be the same and the new
3722  // direction is given by :
3723  if (neighb_pt->Root_pt != Root_pt)
3724  {
3725  int new_up = orig_root_pt->up_equivalent(neighb_pt->Root_pt);
3726 
3727  int new_right = orig_root_pt->right_equivalent(neighb_pt->Root_pt);
3728 
3729  new_dir = rotate(new_up, new_right, direction);
3730  }
3731 
3732  // What's the direction of the interfacial edge when viewed from within
3733  // the neighbour element (including rotations!)
3734  edge = Reflect_edge[new_dir];
3735 
3736  // Get ready to rotate the local coordinates
3737  Vector<double> s_lo_new(3), s_hi_new(3);
3738 
3739  // If the root of the present element is different from the root
3740  // of his neighbour, we have to rotate the lo and hi coordinates
3741  // to have their equivalents from the neighbour's point of view.
3742  if ((neighb_pt->Root_pt != Root_pt))
3743  {
3744  int tmp_dir;
3745  Vector<int> vect1(3);
3746  Vector<int> vect2(3);
3747  Vector<int> vect3(3);
3748  DenseMatrix<int> Mat_rot(3);
3749 
3750  // All this is just to compute the rotation matrix
3751  tmp_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3752  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3753  R);
3754  vect1 = Direction_to_vector[tmp_dir];
3755 
3756 
3757  tmp_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3758  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3759  U);
3760  vect2 = Direction_to_vector[tmp_dir];
3761 
3762 
3763  tmp_dir = rotate(orig_root_pt->up_equivalent(neighb_pt->Root_pt),
3764  orig_root_pt->right_equivalent(neighb_pt->Root_pt),
3765  F);
3766  vect3 = Direction_to_vector[tmp_dir];
3767 
3768 
3769  // Setup the inverse rotation matrix
3770  for (int i = 0; i < 3; i++)
3771  {
3772  Mat_rot(i, 0) = vect1[i];
3773  Mat_rot(i, 1) = vect2[i];
3774  Mat_rot(i, 2) = vect3[i];
3775  }
3776 
3777  // Initialise the translation scheme
3778  Vector<int> translate_s_new(3);
3779 
3780  // Then the rotation of the coordinates
3781  for (unsigned i = 0; i < 3; i++)
3782  {
3783  s_hi_new[i] = 0.0;
3784  s_lo_new[i] = 0.0;
3785  translate_s_new[i] = 0;
3786  for (unsigned k = 0; k < 3; k++)
3787  {
3788  s_hi_new[i] += Mat_rot(i, k) * s_hi[k];
3789  s_lo_new[i] += Mat_rot(i, k) * s_lo[k];
3790  translate_s_new[i] += Mat_rot(i, k) * translate_s[k];
3791  }
3792  }
3793 
3794  s_lo = s_lo_new;
3795  s_hi = s_hi_new;
3796 
3797  // Set the translation scheme
3798  for (unsigned i = 0; i < 3; i++)
3799  {
3800  // abs is ok here; not fabs!
3801  translate_s[i] = std::abs(translate_s_new[i]);
3802  }
3803  }
3804 
3805  } // end if for (neighb_pt!=0)
3806 
3807  return return_pt;
3808  }
oomph::OcTree::edge_neighbour_is_face_neighbour
bool edge_neighbour_is_face_neighbour(const int &edge, OcTree *edge_neighb_pt) const
Definition: octree.cc:2769
oomph::OcTree::Reflect_edge
static Vector< int > Reflect_edge
Get opposite edge, e.g. Reflect_edge[DB]=UF.
Definition: octree.h:335

References abs(), oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, Direct_string, Direction_to_vector, edge_neighbour_is_face_neighbour(), oomph::OcTreeNames::F, gteq_edge_neighbour(), i, k, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::Tree::Level, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, Reflect_edge, oomph::OcTreeNames::RF, oomph::OcTreeRoot::right_equivalent(), oomph::Tree::Root_pt, rotate(), oomph::OcTreeNames::RU, S_base_edge, S_step_edge, RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, oomph::OcTreeNames::UF, and oomph::OcTreeRoot::up_equivalent().

Referenced by doc_true_edge_neighbours(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::node_created_by_neighbour(), oomph::RefineableQElement< 3 >::node_created_by_neighbour(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::node_created_by_son_of_neighbour(), and oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::oc_hang_helper().

◆ mult_mat_mat()

void oomph::OcTree::mult_mat_mat ( const DenseMatrix< int > &  mat1,
const DenseMatrix< int > &  mat2,
DenseMatrix< int > &  mat3 
)
staticprivate

Helper function: Performs the operation : mat3=mat1*mat2.

Helper: Performs the operation Mat3=Mat1*Mat2.

667  {
668  int Sum, i, j, k;
669  for (i = 0; i < 3; i++)
670  {
671  for (j = 0; j < 3; j++)
672  {
673  Sum = 0;
674  for (k = 0; k < 3; k++)
675  {
676  Sum += mat1(i, k) * mat2(k, j);
677  }
678  mat3(i, j) = Sum;
679  }
680  }
681  }
mat1
MatrixXd mat1(size, size)

References i, j, k, and mat1().

Referenced by rotate().

◆ mult_mat_vect()

void oomph::OcTree::mult_mat_vect ( const DenseMatrix< int > &  mat,
const Vector< int > &  vect1,
Vector< int > &  vect2 
)
staticprivate

Helper function: Performs the operation : vect2 = mat*vect1.

Helper: Performs the operation Vect2=Mat*Vect1.

690  {
691  int i, k, sum;
692  for (i = 0; i < 3; i++)
693  {
694  sum = 0;
695  for (k = 0; k < 3; k++)
696  {
697  sum += mat(i, k) * vect1[k];
698  }
699  vect2[i] = sum;
700  }
701  }

References i, and k.

Referenced by rotate().

◆ node_number_to_vertex()

int oomph::OcTree::node_number_to_vertex ( const unsigned &  n,
const unsigned &  nnode1d 
)
static

Return the vertex [LDB,RDB,...] of local (vertex) node n in an element with nnode1d nodes in each coordinate direction.

Return the vertex of local (vertex) node n in an element with nnode1d nodes in each coordinate direction.

492  {
493  using namespace OcTreeNames;
494 
495 #ifdef PARANOID
496  if ((n != 0) && (n != nnode1d - 1) && (n != (nnode1d - 1) * nnode1d) &&
497  (n != nnode1d * nnode1d - 1) &&
498  (n != (nnode1d * nnode1d) * (nnode1d - 1) + 0) &&
499  (n != (nnode1d * nnode1d) * (nnode1d - 1) + nnode1d - 1) &&
500  (n != (nnode1d * nnode1d) * (nnode1d - 1) + (nnode1d - 1) * nnode1d) &&
501  (n != (nnode1d * nnode1d) * (nnode1d - 1) + nnode1d * nnode1d - 1))
502  {
503  std::ostringstream error_stream;
504  error_stream << "Node " << n
505  << " is not a vertex node in a brick element with "
506  << nnode1d << " nodes along each edge!";
507 
508  throw OomphLibError(
509  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
510  }
511 #endif
512 
513  if (n == 0)
514  {
515  return LDB;
516  }
517  else if (n == nnode1d - 1)
518  {
519  return RDB;
520  }
521  else if (n == nnode1d * (nnode1d - 1))
522  {
523  return LUB;
524  }
525  else if (n == nnode1d * nnode1d - 1)
526  {
527  return RUB;
528  }
529  else if (n == nnode1d * nnode1d * (nnode1d - 1))
530  {
531  return LDF;
532  }
533  else if (n == (nnode1d * nnode1d + 1) * (nnode1d - 1))
534  {
535  return RDF;
536  }
537  else if (n == nnode1d * nnode1d * nnode1d - nnode1d)
538  {
539  return LUF;
540  }
541  else if (n == nnode1d * nnode1d * nnode1d - 1)
542  {
543  return RUF;
544  }
545  else
546  {
547  std::ostringstream error_stream;
548  error_stream << "Never get here. local node number: " << n
549  << " is not a vertex node in a brick element with "
550  << nnode1d << " nodes along each edge!" << std::endl;
551  throw OomphLibError(
552  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
553  }
554  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::OcTreeNames::RDF
@ RDF
Definition: octree.h:54
oomph::OcTreeNames::RUB
@ RUB
Definition: octree.h:52
oomph::OcTreeNames::LUF
@ LUF
Definition: octree.h:55
oomph::OcTreeNames::LDF
@ LDF
Definition: octree.h:53
oomph::OcTreeNames::RDB
@ RDB
Definition: octree.h:50
oomph::OcTreeNames::LUB
@ LUB
Definition: octree.h:51
oomph::OcTreeNames::RUF
@ RUF
Definition: octree.h:56
oomph::OcTreeNames::LDB
@ LDB
Definition: octree.h:49

References oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, n, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RUB, and oomph::OcTreeNames::RUF.

Referenced by oomph::OcTreeForest::construct_up_right_equivalents().

◆ operator=()

void oomph::OcTree::operator= ( const OcTree &  )
delete

Broken assignment operator.

◆ rotate() [1/2]

int oomph::OcTree::rotate ( const int &  new_up,
const int &  new_right,
const int &  dir 
)
static

If U[p] becomes new_up and R[ight] becomes new_right then the direction dir becomes rotate(new_up, new_right, dir)

A rotation is defined by the newUp and newRight directions; so if Up becomes newUp and Right becomes newRight then dir becomes rotate(newUp,newRight,dir);

710  {
711  using namespace OcTreeNames;
712 
713 #ifdef PARANOID
714  if ((new_up != L) && (new_up != R) && (new_up != F) && (new_up != B) &&
715  (new_up != U) && (new_up != D))
716  {
717  std::ostringstream error_stream;
718  error_stream << "Wrong new_up: " << Direct_string[new_up] << std::endl;
719  throw OomphLibError(
720  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
721  }
722  if ((new_right != L) && (new_right != R) && (new_right != F) &&
723  (new_right != B) && (new_right != U) && (new_right != D))
724  {
725  std::ostringstream error_stream;
726  error_stream << "Wrong new_right: " << Direct_string[new_right]
727  << std::endl;
728  throw OomphLibError(
729  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
730  }
731 #endif
732 
733  Vector<int> vect_dir(3);
734  Vector<int> vect_new_dir(3);
735  // translate the direction to a vector
736  vect_dir = Direction_to_vector[dir];
737 
738  // rotate it
739  vect_new_dir = rotate(new_up, new_right, vect_dir);
740 
741  // translate the vector back to a direction
742  return Vector_to_direction[vect_new_dir];
743  }

References D, Direct_string, Direction_to_vector, oomph::OcTreeNames::F, L, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, rotate(), RachelsAdvectionDiffusion::U, and Vector_to_direction.

◆ rotate() [2/2]

Vector< int > oomph::OcTree::rotate ( const int &  new_up,
const int &  new_right,
const Vector< int > &  dir 
)
static

If U[p] becomes new_up and R[ight] becomes new_right then the direction vector dir becomes rotate(new_up, new_right, dir)

This function rotates a vector according to a rotation of the axes that changes up to new_up and right to new_right.

753  {
754  using namespace OcTreeNames;
755 
756 #ifdef PARANOID
757  if ((new_up != L) && (new_up != R) && (new_up != F) && (new_up != B) &&
758  (new_up != U) && (new_up != D))
759  {
760  std::ostringstream error_stream;
761  error_stream << "Wrong new_up: " << Direct_string[new_up] << std::endl;
762  throw OomphLibError(
763  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
764  }
765  if ((new_right != L) && (new_right != R) && (new_right != F) &&
766  (new_right != B) && (new_right != U) && (new_right != D))
767  {
768  std::ostringstream error_stream;
769  error_stream << "Wrong new_right: " << Direct_string[new_right]
770  << std::endl;
771  throw OomphLibError(
772  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
773  }
774 #endif
775 
776  // Every possible rotation of an element can be produced by the
777  // compostition of two (or one) rotations about one of the main
778  // axes of the element (R, U, F). So for each (newRight,newUp)
779  // we define the (angle1,axis1) of the first rotation and the
780  // (angle2, axis2) of the second one (if needed)
781 
782  int axis1, axis2, angle1, angle2, nrot;
783  nrot = 2;
784  if (new_up == U)
785  {
786  nrot = 1;
787  axis1 = U;
788  switch (new_right)
789  {
790  case R:
791  angle1 = 0;
792  break;
793  case B:
794  angle1 = 1;
795  break;
796  case L:
797  angle1 = 2;
798  break;
799  case F:
800  angle1 = 3;
801  break;
802  default:
803  std::ostringstream error_stream;
804  error_stream << "New_right is " << new_right << " ("
805  << Direct_string[new_right] << "). "
806  << "It should be R, B, L, or F" << std::endl;
807  throw OomphLibError(error_stream.str(),
808  OOMPH_CURRENT_FUNCTION,
809  OOMPH_EXCEPTION_LOCATION);
810  }
811  }
812 
813  if (new_up == D)
814  {
815  switch (new_right)
816  {
817  case R:
818  nrot = 1;
819  axis1 = R;
820  angle1 = 2;
821  break;
822  case B:
823  axis1 = R;
824  angle1 = 2;
825  axis2 = U;
826  angle2 = 1;
827  break;
828  case L:
829  nrot = 1;
830  axis1 = F;
831  angle1 = 2;
832  break;
833  case F:
834  axis1 = R;
835  angle1 = 2;
836  axis2 = U;
837  angle2 = 3;
838  break;
839  default:
840  std::ostringstream error_stream;
841  error_stream << "New_right is " << new_right << " ("
842  << Direct_string[new_right] << "). "
843  << "It should be R, B, L, or F" << std::endl;
844  throw OomphLibError(error_stream.str(),
845  OOMPH_CURRENT_FUNCTION,
846  OOMPH_EXCEPTION_LOCATION);
847  }
848  }
849 
850  if (new_up == R)
851  {
852  switch (new_right)
853  {
854  case D:
855  nrot = 1;
856  axis1 = F;
857  angle1 = 3;
858  break;
859  case B:
860  axis1 = F;
861  angle1 = 3;
862  axis2 = R;
863  angle2 = 1;
864  break;
865  case U:
866  axis1 = F;
867  angle1 = 1;
868  axis2 = U;
869  angle2 = 2;
870  break;
871  case F:
872  axis1 = F;
873  angle1 = 3;
874  axis2 = R;
875  angle2 = 3;
876  break;
877  default:
878  std::ostringstream error_stream;
879  error_stream << "New_right is " << new_right << " ("
880  << Direct_string[new_right] << "). "
881  << "It should be D, B, U, or F" << std::endl;
882  throw OomphLibError(error_stream.str(),
883  OOMPH_CURRENT_FUNCTION,
884  OOMPH_EXCEPTION_LOCATION);
885  }
886  }
887 
888  if (new_up == L)
889  {
890  switch (new_right)
891  {
892  case D:
893  axis1 = F;
894  angle1 = 1;
895  axis2 = R;
896  angle2 = 2;
897  break;
898  case B:
899  axis1 = F;
900  angle1 = 1;
901  axis2 = R;
902  angle2 = 3;
903  break;
904  case U:
905  nrot = 1;
906  axis1 = F;
907  angle1 = 1;
908  break;
909  case F:
910  axis1 = F;
911  angle1 = 1;
912  axis2 = R;
913  angle2 = 1;
914  break;
915  default:
916  std::ostringstream error_stream;
917  error_stream << "New_right is " << new_right << " ("
918  << Direct_string[new_right] << "). "
919  << "It should be D, B, U, or F" << std::endl;
920  throw OomphLibError(error_stream.str(),
921  OOMPH_CURRENT_FUNCTION,
922  OOMPH_EXCEPTION_LOCATION);
923  }
924  }
925 
926  if (new_up == F)
927  {
928  switch (new_right)
929  {
930  case R:
931  nrot = 1;
932  axis1 = R;
933  angle1 = 1;
934  break;
935  case L:
936  axis1 = R;
937  angle1 = 1;
938  axis2 = F;
939  angle2 = 2;
940  break;
941  case U:
942  axis1 = R;
943  angle1 = 1;
944  axis2 = F;
945  angle2 = 1;
946  break;
947  case D:
948  axis1 = R;
949  angle1 = 1;
950  axis2 = F;
951  angle2 = 3;
952  break;
953  default:
954  std::ostringstream error_stream;
955  error_stream << "New_right is " << new_right << " ("
956  << Direct_string[new_right] << "). "
957  << "It should be R, L, U, or D" << std::endl;
958  throw OomphLibError(error_stream.str(),
959  OOMPH_CURRENT_FUNCTION,
960  OOMPH_EXCEPTION_LOCATION);
961  }
962  }
963  if (new_up == B)
964  {
965  switch (new_right)
966  {
967  case R:
968  nrot = 1;
969  axis1 = R;
970  angle1 = 3;
971  break;
972  case L:
973  axis1 = R;
974  angle1 = 3;
975  axis2 = F;
976  angle2 = 2;
977  break;
978  case U:
979  axis1 = R;
980  angle1 = 3;
981  axis2 = F;
982  angle2 = 1;
983  break;
984  case D:
985  axis1 = R;
986  angle1 = 3;
987  axis2 = F;
988  angle2 = 3;
989  break;
990  default:
991  std::ostringstream error_stream;
992  error_stream << "New_right is " << new_right << " ("
993  << Direct_string[new_right] << "). "
994  << "It should be R, L, U, or D" << std::endl;
995  throw OomphLibError(error_stream.str(),
996  OOMPH_CURRENT_FUNCTION,
997  OOMPH_EXCEPTION_LOCATION);
998  }
999  }
1000 
1001  Vector<int> vect_new_dir(3);
1002  DenseMatrix<int> mat1(3);
1003  DenseMatrix<int> mat2(3);
1004  DenseMatrix<int> mat3(3);
1005 
1006  // Then we build the first rotation matrix
1007  construct_rotation_matrix(axis1, angle1, mat1);
1008 
1009  // If needed we build the second one
1010  if (nrot == 2)
1011  {
1012  // And we make the composition of the two rotations
1013  // which is stored in Mat3
1014  construct_rotation_matrix(axis2, angle2, mat2);
1015  mult_mat_mat(mat2, mat1, mat3);
1016  }
1017  else
1018  {
1019  // Else we just copy Mat1 into Mat3
1020  for (int i = 0; i < 3; i++)
1021  {
1022  for (int j = 0; j < 3; j++)
1023  {
1024  mat3(i, j) = mat1(i, j);
1025  }
1026  }
1027  }
1028 
1029  // Rotate
1030  mult_mat_vect(mat3, dir, vect_new_dir);
1031 
1032  // Return the corresponding vector
1033  return vect_new_dir;
1034  }
oomph::OcTree::mult_mat_vect
static void mult_mat_vect(const DenseMatrix< int > &mat, const Vector< int > &vect1, Vector< int > &vect2)
Helper function: Performs the operation : vect2 = mat*vect1.
Definition: octree.cc:687
oomph::OcTree::construct_rotation_matrix
static void construct_rotation_matrix(int &axis, int &angle, DenseMatrix< int > &mat)
Definition: octree.cc:609
oomph::OcTree::mult_mat_mat
static void mult_mat_mat(const DenseMatrix< int > &mat1, const DenseMatrix< int > &mat2, DenseMatrix< int > &mat3)
Helper function: Performs the operation : mat3=mat1*mat2.
Definition: octree.cc:664

References construct_rotation_matrix(), D, Direct_string, oomph::OcTreeNames::F, i, j, L, mat1(), mult_mat_mat(), mult_mat_vect(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, and RachelsAdvectionDiffusion::U.

Referenced by gteq_edge_neighbour(), gteq_face_neighbour(), gteq_true_edge_neighbour(), rotate(), and setup_static_data().

◆ self_test()

unsigned oomph::OcTree::self_test ( )

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.

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 neigbour. . If the difference is less than Tree::Max_neighbour_finding_tolerance. return success (0), otherwise failure (1)

4201  {
4202  // Stick pointers to all nodes into Vector and number elements
4203  // in the process
4204  Vector<Tree*> all_nodes_pt;
4205  stick_all_tree_nodes_into_vector(all_nodes_pt);
4206 
4207  long int count = 0;
4208  unsigned long num_nodes = all_nodes_pt.size();
4209 
4210  for (unsigned long i = 0; i < num_nodes; i++)
4211  {
4212  all_nodes_pt[i]->object_pt()->set_number(++count);
4213  }
4214 
4215  // Check neighbours vertices and their opposite points
4216  std::ofstream neighbours_file;
4217  std::ofstream no_true_edge_file;
4218  std::ofstream neighbours_txt_file;
4219 
4220  double max_error_face = 0.0;
4221  OcTree::doc_face_neighbours(
4222  all_nodes_pt, neighbours_file, neighbours_txt_file, max_error_face);
4223 
4224  double max_error_edge = 0.0;
4225  OcTree::doc_true_edge_neighbours(all_nodes_pt,
4226  neighbours_file,
4227  no_true_edge_file,
4228  neighbours_txt_file,
4229  max_error_edge);
4230  bool failed = false;
4231  if (max_error_face > Max_neighbour_finding_tolerance)
4232  {
4233  oomph_info
4234  << "\n \n Failed self_test() for OcTree because of faces: Max. error "
4235  << max_error_face << std::endl
4236  << std::endl;
4237  failed = true;
4238  }
4239 
4240  if (max_error_edge > Max_neighbour_finding_tolerance)
4241  {
4242  oomph_info
4243  << "\n \n Failed self_test() for OcTree because of edges: Max. error "
4244  << max_error_edge << std::endl
4245  << std::endl;
4246  failed = true;
4247  }
4248 
4249  double max_error = max_error_face;
4250  if (max_error_edge > max_error) max_error = max_error_edge;
4251 
4252  if (failed)
4253  {
4254  return 1;
4255  }
4256  else
4257  {
4258  oomph_info << "Passed self_test() for OcTree: Max. error " << max_error
4259  << std::endl;
4260  return 0;
4261  }
4262  }
oomph::OcTree::doc_true_edge_neighbours
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)
Definition: octree.cc:4552
oomph::OcTree::doc_face_neighbours
static void doc_face_neighbours(Vector< Tree * > forest_nodes_pt, std::ofstream &neighbours_file, std::ofstream &neighbours_txt_file, double &max_error)
Definition: octree.cc:4275
oomph::Tree::stick_all_tree_nodes_into_vector
void stick_all_tree_nodes_into_vector(Vector< Tree * > &)
Traverse and stick pointers to all "nodes" into Vector.
Definition: tree.cc:277
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References doc_face_neighbours(), doc_true_edge_neighbours(), i, MeshRefinement::max_error, oomph::Tree::Max_neighbour_finding_tolerance, oomph::oomph_info, and oomph::Tree::stick_all_tree_nodes_into_vector().

◆ setup_static_data()

void oomph::OcTree::setup_static_data ( )
static

Setup the static data, rotation and reflection schemes, etc.

Setup static data for OcTree.

1041  {
1042  using namespace OcTreeNames;
1043 
1044 
1045 #ifdef PARANOID
1046  if (Tree::OMEGA != OcTree::OMEGA)
1047  {
1048  std::ostringstream error_stream;
1049  error_stream << "Inconsistent enumeration! \n Tree::OMEGA="
1050  << Tree::OMEGA << "\nOcTree::OMEGA=" << OcTree::OMEGA
1051  << std::endl;
1052  throw OomphLibError(
1053  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1054  }
1055 #endif
1056 
1057  // Private static data
1058  //--------------------
1059 
1060  // Set flag to indicate that static data has been setup
1061  Static_data_has_been_setup = true;
1062 
1063  // Initialisation of cosi and sini used in rotation matrix
1064  Cosi.resize(4);
1065  Sini.resize(4);
1066 
1067  Cosi[0] = 1;
1068  Cosi[1] = 0;
1069  Cosi[2] = -1;
1070  Cosi[3] = 0;
1071 
1072  Sini[0] = 0;
1073  Sini[1] = 1;
1074  Sini[2] = 0;
1075  Sini[3] = -1;
1076 
1077 
1078  // Build direction/octant adjacency scheme
1079  // Is_adjacent(i_face,j_octant):
1080  // Is face adjacent to octant?
1081  // See the table in Samet book
1082  Is_adjacent.resize(27, 27);
1083  Is_adjacent(L, LDB) = true;
1084  Is_adjacent(R, LDB) = false;
1085  Is_adjacent(D, LDB) = true;
1086  Is_adjacent(U, LDB) = false;
1087  Is_adjacent(B, LDB) = true;
1088  Is_adjacent(F, LDB) = false;
1089 
1090  Is_adjacent(L, LDF) = true;
1091  Is_adjacent(R, LDF) = false;
1092  Is_adjacent(D, LDF) = true;
1093  Is_adjacent(U, LDF) = false;
1094  Is_adjacent(B, LDF) = false;
1095  Is_adjacent(F, LDF) = true;
1096 
1097  Is_adjacent(L, LUB) = true;
1098  Is_adjacent(R, LUB) = false;
1099  Is_adjacent(D, LUB) = false;
1100  Is_adjacent(U, LUB) = true;
1101  Is_adjacent(B, LUB) = true;
1102  Is_adjacent(F, LUB) = false;
1103 
1104  Is_adjacent(L, LUF) = true;
1105  Is_adjacent(R, LUF) = false;
1106  Is_adjacent(D, LUF) = false;
1107  Is_adjacent(U, LUF) = true;
1108  Is_adjacent(B, LUF) = false;
1109  Is_adjacent(F, LUF) = true;
1110 
1111  Is_adjacent(L, RDB) = false;
1112  Is_adjacent(R, RDB) = true;
1113  Is_adjacent(D, RDB) = true;
1114  Is_adjacent(U, RDB) = false;
1115  Is_adjacent(B, RDB) = true;
1116  Is_adjacent(F, RDB) = false;
1117 
1118  Is_adjacent(L, RDF) = false;
1119  Is_adjacent(R, RDF) = true;
1120  Is_adjacent(D, RDF) = true;
1121  Is_adjacent(U, RDF) = false;
1122  Is_adjacent(B, RDF) = false;
1123  Is_adjacent(F, RDF) = true;
1124 
1125  Is_adjacent(L, RUB) = false;
1126  Is_adjacent(R, RUB) = true;
1127  Is_adjacent(D, RUB) = false;
1128  Is_adjacent(U, RUB) = true;
1129  Is_adjacent(B, RUB) = true;
1130  Is_adjacent(F, RUB) = false;
1131 
1132  Is_adjacent(L, RUF) = false;
1133  Is_adjacent(R, RUF) = true;
1134  Is_adjacent(D, RUF) = false;
1135  Is_adjacent(U, RUF) = true;
1136  Is_adjacent(B, RUF) = false;
1137  Is_adjacent(F, RUF) = true;
1138 
1139 
1140  // Build direction/octant adjacency scheme
1141  // Is_adjacent(i_edge,j_octant):
1142  // Is edge adjacent to octant?
1143  // See the table in Samet book
1144  Is_adjacent(LD, LDB) = true;
1145  Is_adjacent(LD, LDF) = true;
1146  Is_adjacent(LD, LUB) = false;
1147  Is_adjacent(LD, LUF) = false;
1148  Is_adjacent(LD, RDB) = false;
1149  Is_adjacent(LD, RDF) = false;
1150  Is_adjacent(LD, RUB) = false;
1151  Is_adjacent(LD, RUF) = false;
1152 
1153  Is_adjacent(LU, LDB) = false;
1154  Is_adjacent(LU, LDF) = false;
1155  Is_adjacent(LU, LUB) = true;
1156  Is_adjacent(LU, LUF) = true;
1157  Is_adjacent(LU, RDB) = false;
1158  Is_adjacent(LU, RDF) = false;
1159  Is_adjacent(LU, RUB) = false;
1160  Is_adjacent(LU, RUF) = false;
1161 
1162 
1163  Is_adjacent(LB, LDB) = true;
1164  Is_adjacent(LB, LDF) = false;
1165  Is_adjacent(LB, LUB) = true;
1166  Is_adjacent(LB, LUF) = false;
1167  Is_adjacent(LB, RDB) = false;
1168  Is_adjacent(LB, RDF) = false;
1169  Is_adjacent(LB, RUB) = false;
1170  Is_adjacent(LB, RUF) = false;
1171 
1172 
1173  Is_adjacent(LF, LDB) = false;
1174  Is_adjacent(LF, LDF) = true;
1175  Is_adjacent(LF, LUB) = false;
1176  Is_adjacent(LF, LUF) = true;
1177  Is_adjacent(LF, RDB) = false;
1178  Is_adjacent(LF, RDF) = false;
1179  Is_adjacent(LF, RUB) = false;
1180  Is_adjacent(LF, RUF) = false;
1181 
1182 
1183  Is_adjacent(RD, LDB) = false;
1184  Is_adjacent(RD, LDF) = false;
1185  Is_adjacent(RD, LUB) = false;
1186  Is_adjacent(RD, LUF) = false;
1187  Is_adjacent(RD, RDB) = true;
1188  Is_adjacent(RD, RDF) = true;
1189  Is_adjacent(RD, RUB) = false;
1190  Is_adjacent(RD, RUF) = false;
1191 
1192 
1193  Is_adjacent(RU, LDB) = false;
1194  Is_adjacent(RU, LDF) = false;
1195  Is_adjacent(RU, LUB) = false;
1196  Is_adjacent(RU, LUF) = false;
1197  Is_adjacent(RU, RDB) = false;
1198  Is_adjacent(RU, RDF) = false;
1199  Is_adjacent(RU, RUB) = true;
1200  Is_adjacent(RU, RUF) = true;
1201 
1202  Is_adjacent(RB, LDB) = false;
1203  Is_adjacent(RB, LDF) = false;
1204  Is_adjacent(RB, LUB) = false;
1205  Is_adjacent(RB, LUF) = false;
1206  Is_adjacent(RB, RDB) = true;
1207  Is_adjacent(RB, RDF) = false;
1208  Is_adjacent(RB, RUB) = true;
1209  Is_adjacent(RB, RUF) = false;
1210 
1211  Is_adjacent(RF, LDB) = false;
1212  Is_adjacent(RF, LDF) = false;
1213  Is_adjacent(RF, LUB) = false;
1214  Is_adjacent(RF, LUF) = false;
1215  Is_adjacent(RF, RDB) = false;
1216  Is_adjacent(RF, RDF) = true;
1217  Is_adjacent(RF, RUB) = false;
1218  Is_adjacent(RF, RUF) = true;
1219 
1220 
1221  Is_adjacent(DB, LDB) = true;
1222  Is_adjacent(DB, LDF) = false;
1223  Is_adjacent(DB, LUB) = false;
1224  Is_adjacent(DB, LUF) = false;
1225  Is_adjacent(DB, RDB) = true;
1226  Is_adjacent(DB, RDF) = false;
1227  Is_adjacent(DB, RUB) = false;
1228  Is_adjacent(DB, RUF) = false;
1229 
1230 
1231  Is_adjacent(DF, LDB) = false;
1232  Is_adjacent(DF, LDF) = true;
1233  Is_adjacent(DF, LUB) = false;
1234  Is_adjacent(DF, LUF) = false;
1235  Is_adjacent(DF, RDB) = false;
1236  Is_adjacent(DF, RDF) = true;
1237  Is_adjacent(DF, RUB) = false;
1238  Is_adjacent(DF, RUF) = false;
1239 
1240  Is_adjacent(UB, LDB) = false;
1241  Is_adjacent(UB, LDF) = false;
1242  Is_adjacent(UB, LUB) = true;
1243  Is_adjacent(UB, LUF) = false;
1244  Is_adjacent(UB, RDB) = false;
1245  Is_adjacent(UB, RDF) = false;
1246  Is_adjacent(UB, RUB) = true;
1247  Is_adjacent(UB, RUF) = false;
1248 
1249  Is_adjacent(UF, LDB) = false;
1250  Is_adjacent(UF, LDF) = false;
1251  Is_adjacent(UF, LUB) = false;
1252  Is_adjacent(UF, LUF) = true;
1253  Is_adjacent(UF, RDB) = false;
1254  Is_adjacent(UF, RDF) = false;
1255  Is_adjacent(UF, RUB) = false;
1256  Is_adjacent(UF, RUF) = true;
1257 
1258 
1259  // Common face of various octants from Samet's book
1260  Common_face.resize(27, 27);
1261  Common_face(LDB, LDB) = OMEGA;
1262  Common_face(LDB, LDF) = OMEGA;
1263  Common_face(LDB, LUB) = OMEGA;
1264  Common_face(LDB, LUF) = L;
1265  Common_face(LDB, RDB) = OMEGA;
1266  Common_face(LDB, RDF) = D;
1267  Common_face(LDB, RUB) = B;
1268  Common_face(LDB, RUF) = OMEGA;
1269 
1270  Common_face(LDF, LDB) = OMEGA;
1271  Common_face(LDF, LDF) = OMEGA;
1272  Common_face(LDF, LUB) = L;
1273  Common_face(LDF, LUF) = OMEGA;
1274  Common_face(LDF, RDB) = D;
1275  Common_face(LDF, RDF) = OMEGA;
1276  Common_face(LDF, RUB) = OMEGA;
1277  Common_face(LDF, RUF) = F;
1278 
1279  Common_face(LUB, LDB) = OMEGA;
1280  Common_face(LUB, LDF) = L;
1281  Common_face(LUB, LUB) = OMEGA;
1282  Common_face(LUB, LUF) = OMEGA;
1283  Common_face(LUB, RDB) = B;
1284  Common_face(LUB, RDF) = OMEGA;
1285  Common_face(LUB, RUB) = OMEGA;
1286  Common_face(LUB, RUF) = U;
1287 
1288  Common_face(LUF, LDB) = L;
1289  Common_face(LUF, LDF) = OMEGA;
1290  Common_face(LUF, LUB) = OMEGA;
1291  Common_face(LUF, LUF) = OMEGA;
1292  Common_face(LUF, RDB) = OMEGA;
1293  Common_face(LUF, RDF) = F;
1294  Common_face(LUF, RUB) = U;
1295  Common_face(LUF, RUF) = OMEGA;
1296 
1297  Common_face(RDB, LDB) = OMEGA;
1298  Common_face(RDB, LDF) = D;
1299  Common_face(RDB, LUB) = B;
1300  Common_face(RDB, LUF) = OMEGA;
1301  Common_face(RDB, RDB) = OMEGA;
1302  Common_face(RDB, RDF) = OMEGA;
1303  Common_face(RDB, RUB) = OMEGA;
1304  Common_face(RDB, RUF) = R;
1305 
1306  Common_face(RDF, LDB) = D;
1307  Common_face(RDF, LDF) = OMEGA;
1308  Common_face(RDF, LUB) = OMEGA;
1309  Common_face(RDF, LUF) = F;
1310  Common_face(RDF, RDB) = OMEGA;
1311  Common_face(RDF, RDF) = OMEGA;
1312  Common_face(RDF, RUB) = R;
1313  Common_face(RDF, RUF) = OMEGA;
1314 
1315  Common_face(RUB, LDB) = B;
1316  Common_face(RUB, LDF) = OMEGA;
1317  Common_face(RUB, LUB) = OMEGA;
1318  Common_face(RUB, LUF) = U;
1319  Common_face(RUB, RDB) = OMEGA;
1320  Common_face(RUB, RDF) = R;
1321  Common_face(RUB, RUB) = OMEGA;
1322  Common_face(RUB, RUF) = OMEGA;
1323 
1324  Common_face(RUF, LDB) = OMEGA;
1325  Common_face(RUF, LDF) = F;
1326  Common_face(RUF, LUB) = U;
1327  Common_face(RUF, LUF) = OMEGA;
1328  Common_face(RUF, RDB) = R;
1329  Common_face(RUF, RDF) = OMEGA;
1330  Common_face(RUF, RUB) = OMEGA;
1331  Common_face(RUF, RUF) = OMEGA;
1332 
1333 
1334  // Common face of various edges/octants from Samet's book
1335  Common_face(LD, LDB) = OMEGA;
1336  Common_face(LD, LDF) = OMEGA;
1337  Common_face(LD, LUB) = L;
1338  Common_face(LD, LUF) = L;
1339  Common_face(LD, RDB) = D;
1340  Common_face(LD, RDF) = D;
1341  Common_face(LD, RUB) = OMEGA;
1342  Common_face(LD, RUF) = OMEGA;
1343 
1344  Common_face(LU, LDB) = L;
1345  Common_face(LU, LDF) = L;
1346  Common_face(LU, LUB) = OMEGA;
1347  Common_face(LU, LUF) = OMEGA;
1348  Common_face(LU, RDB) = OMEGA;
1349  Common_face(LU, RDF) = OMEGA;
1350  Common_face(LU, RUB) = U;
1351  Common_face(LU, RUF) = U;
1352 
1353  Common_face(LB, LDB) = OMEGA;
1354  Common_face(LB, LDF) = L;
1355  Common_face(LB, LUB) = OMEGA;
1356  Common_face(LB, LUF) = L;
1357  Common_face(LB, RDB) = B;
1358  Common_face(LB, RDF) = OMEGA;
1359  Common_face(LB, RUB) = B;
1360  Common_face(LB, RUF) = OMEGA;
1361 
1362  Common_face(LF, LDB) = L;
1363  Common_face(LF, LDF) = OMEGA;
1364  Common_face(LF, LUB) = L;
1365  Common_face(LF, LUF) = OMEGA;
1366  Common_face(LF, RDB) = OMEGA;
1367  Common_face(LF, RDF) = F;
1368  Common_face(LF, RUB) = OMEGA;
1369  Common_face(LF, RUF) = F;
1370 
1371  Common_face(RD, LDB) = D;
1372  Common_face(RD, LDF) = D;
1373  Common_face(RD, LUB) = OMEGA;
1374  Common_face(RD, LUF) = OMEGA;
1375  Common_face(RD, RDB) = OMEGA;
1376  Common_face(RD, RDF) = OMEGA;
1377  Common_face(RD, RUB) = R;
1378  Common_face(RD, RUF) = R;
1379 
1380  Common_face(RU, LDB) = OMEGA;
1381  Common_face(RU, LDF) = OMEGA;
1382  Common_face(RU, LUB) = U;
1383  Common_face(RU, LUF) = U;
1384  Common_face(RU, RDB) = R;
1385  Common_face(RU, RDF) = R;
1386  Common_face(RU, RUB) = OMEGA;
1387  Common_face(RU, RUF) = OMEGA;
1388 
1389  Common_face(RB, LDB) = B;
1390  Common_face(RB, LDF) = OMEGA;
1391  Common_face(RB, LUB) = B;
1392  Common_face(RB, LUF) = OMEGA;
1393  Common_face(RB, RDB) = OMEGA;
1394  Common_face(RB, RDF) = R;
1395  Common_face(RB, RUB) = OMEGA;
1396  Common_face(RB, RUF) = R;
1397 
1398  Common_face(RF, LDB) = OMEGA;
1399  Common_face(RF, LDF) = F;
1400  Common_face(RF, LUB) = OMEGA;
1401  Common_face(RF, LUF) = F;
1402  Common_face(RF, RDB) = R;
1403  Common_face(RF, RDF) = OMEGA;
1404  Common_face(RF, RUB) = R;
1405  Common_face(RF, RUF) = OMEGA;
1406 
1407 
1408  Common_face(DB, LDB) = OMEGA;
1409  Common_face(DB, LDF) = D;
1410  Common_face(DB, LUB) = B;
1411  Common_face(DB, LUF) = OMEGA;
1412  Common_face(DB, RDB) = OMEGA;
1413  Common_face(DB, RDF) = D;
1414  Common_face(DB, RUB) = B;
1415  Common_face(DB, RUF) = OMEGA;
1416 
1417  Common_face(DF, LDB) = D;
1418  Common_face(DF, LDF) = OMEGA;
1419  Common_face(DF, LUB) = OMEGA;
1420  Common_face(DF, LUF) = F;
1421  Common_face(DF, RDB) = D;
1422  Common_face(DF, RDF) = OMEGA;
1423  Common_face(DF, RUB) = OMEGA;
1424  Common_face(DF, RUF) = F;
1425 
1426  Common_face(UB, LDB) = B;
1427  Common_face(UB, LDF) = OMEGA;
1428  Common_face(UB, LUB) = OMEGA;
1429  Common_face(UB, LUF) = U;
1430  Common_face(UB, RDB) = B;
1431  Common_face(UB, RDF) = OMEGA;
1432  Common_face(UB, RUB) = OMEGA;
1433  Common_face(UB, RUF) = U;
1434 
1435  Common_face(UF, LDB) = OMEGA;
1436  Common_face(UF, LDF) = F;
1437  Common_face(UF, LUB) = U;
1438  Common_face(UF, LUF) = OMEGA;
1439  Common_face(UF, RDB) = OMEGA;
1440  Common_face(UF, RDF) = F;
1441  Common_face(UF, RUB) = U;
1442  Common_face(UF, RUF) = OMEGA;
1443 
1444 
1445  // Tecplot colours for neighbours in various directions
1446  Colour.resize(27);
1447  Colour[R] = "CYAN";
1448  Colour[L] = "RED";
1449  Colour[U] = "GREEN";
1450  Colour[D] = "BLUE";
1451  Colour[F] = "CUSTOM3";
1452  Colour[B] = "PURPLE";
1453  Colour[OMEGA] = "YELLOW";
1454 
1455  Colour[LB] = "BLUE";
1456  Colour[RB] = "BLUE";
1457  Colour[DB] = "BLUE";
1458  Colour[UB] = "BLUE";
1459 
1460  Colour[LD] = "GREEN";
1461  Colour[RD] = "GREEN";
1462  Colour[LU] = "GREEN";
1463  Colour[RU] = "GREEN";
1464 
1465 
1466  Colour[LF] = "RED";
1467  Colour[RF] = "RED";
1468  Colour[DF] = "RED";
1469  Colour[UF] = "RED";
1470 
1471  Colour[OMEGA] = "YELLOW";
1472 
1473 
1474  // Reflection scheme:
1475  // Reflect(direction,octant): Get mirror of octant in direction
1476  // Table in Samet book as well
1477  Reflect.resize(27, 27);
1478 
1479  Reflect(L, LDB) = RDB;
1480  Reflect(R, LDB) = RDB;
1481  Reflect(D, LDB) = LUB;
1482  Reflect(U, LDB) = LUB;
1483  Reflect(B, LDB) = LDF;
1484  Reflect(F, LDB) = LDF;
1485 
1486  Reflect(L, LDF) = RDF;
1487  Reflect(R, LDF) = RDF;
1488  Reflect(D, LDF) = LUF;
1489  Reflect(U, LDF) = LUF;
1490  Reflect(B, LDF) = LDB;
1491  Reflect(F, LDF) = LDB;
1492 
1493  Reflect(L, LUB) = RUB;
1494  Reflect(R, LUB) = RUB;
1495  Reflect(D, LUB) = LDB;
1496  Reflect(U, LUB) = LDB;
1497  Reflect(B, LUB) = LUF;
1498  Reflect(F, LUB) = LUF;
1499 
1500  Reflect(L, LUF) = RUF;
1501  Reflect(R, LUF) = RUF;
1502  Reflect(D, LUF) = LDF;
1503  Reflect(U, LUF) = LDF;
1504  Reflect(B, LUF) = LUB;
1505  Reflect(F, LUF) = LUB;
1506 
1507  Reflect(L, RDB) = LDB;
1508  Reflect(R, RDB) = LDB;
1509  Reflect(D, RDB) = RUB;
1510  Reflect(U, RDB) = RUB;
1511  Reflect(B, RDB) = RDF;
1512  Reflect(F, RDB) = RDF;
1513 
1514  Reflect(L, RDF) = LDF;
1515  Reflect(R, RDF) = LDF;
1516  Reflect(D, RDF) = RUF;
1517  Reflect(U, RDF) = RUF;
1518  Reflect(B, RDF) = RDB;
1519  Reflect(F, RDF) = RDB;
1520 
1521  Reflect(L, RUB) = LUB;
1522  Reflect(R, RUB) = LUB;
1523  Reflect(D, RUB) = RDB;
1524  Reflect(U, RUB) = RDB;
1525  Reflect(B, RUB) = RUF;
1526  Reflect(F, RUB) = RUF;
1527 
1528  Reflect(L, RUF) = LUF;
1529  Reflect(R, RUF) = LUF;
1530  Reflect(D, RUF) = RDF;
1531  Reflect(U, RUF) = RDF;
1532  Reflect(B, RUF) = RUB;
1533  Reflect(F, RUF) = RUB;
1534 
1535 
1536  // Reflection scheme:
1537  // Reflect(direction (edge) ,octant): Get mirror of edge in direction
1538  // Table in Samet book as well
1539 
1540  Reflect(LD, LDB) = RUB;
1541  Reflect(LU, LDB) = RUB;
1542  Reflect(RD, LDB) = RUB;
1543  Reflect(RU, LDB) = RUB;
1544 
1545  Reflect(LD, LDF) = RUF;
1546  Reflect(LU, LDF) = RUF;
1547  Reflect(RD, LDF) = RUF;
1548  Reflect(RU, LDF) = RUF;
1549 
1550  Reflect(LD, LUB) = RDB;
1551  Reflect(LU, LUB) = RDB;
1552  Reflect(RD, LUB) = RDB;
1553  Reflect(RU, LUB) = RDB;
1554 
1555  Reflect(LD, LUF) = RDF;
1556  Reflect(LU, LUF) = RDF;
1557  Reflect(RD, LUF) = RDF;
1558  Reflect(RU, LUF) = RDF;
1559 
1560 
1561  Reflect(LD, RDB) = LUB;
1562  Reflect(LU, RDB) = LUB;
1563  Reflect(RD, RDB) = LUB;
1564  Reflect(RU, RDB) = LUB;
1565 
1566  Reflect(LD, RDF) = LUF;
1567  Reflect(LU, RDF) = LUF;
1568  Reflect(RD, RDF) = LUF;
1569  Reflect(RU, RDF) = LUF;
1570 
1571  Reflect(LD, RUB) = LDB;
1572  Reflect(LU, RUB) = LDB;
1573  Reflect(RD, RUB) = LDB;
1574  Reflect(RU, RUB) = LDB;
1575 
1576  Reflect(LD, RUF) = LDF;
1577  Reflect(LU, RUF) = LDF;
1578  Reflect(RD, RUF) = LDF;
1579  Reflect(RU, RUF) = LDF;
1580 
1581 
1582  Reflect(LB, LDB) = RDF;
1583  Reflect(LF, LDB) = RDF;
1584  Reflect(RB, LDB) = RDF;
1585  Reflect(RF, LDB) = RDF;
1586 
1587  Reflect(LB, LDF) = RDB;
1588  Reflect(LF, LDF) = RDB;
1589  Reflect(RB, LDF) = RDB;
1590  Reflect(RF, LDF) = RDB;
1591 
1592  Reflect(LB, LUB) = RUF;
1593  Reflect(LF, LUB) = RUF;
1594  Reflect(RB, LUB) = RUF;
1595  Reflect(RF, LUB) = RUF;
1596 
1597  Reflect(LB, LUF) = RUB;
1598  Reflect(LF, LUF) = RUB;
1599  Reflect(RB, LUF) = RUB;
1600  Reflect(RF, LUF) = RUB;
1601 
1602  Reflect(LB, RDB) = LDF;
1603  Reflect(LF, RDB) = LDF;
1604  Reflect(RB, RDB) = LDF;
1605  Reflect(RF, RDB) = LDF;
1606 
1607  Reflect(LB, RDF) = LDB;
1608  Reflect(LF, RDF) = LDB;
1609  Reflect(RB, RDF) = LDB;
1610  Reflect(RF, RDF) = LDB;
1611 
1612  Reflect(LB, RUB) = LUF;
1613  Reflect(LF, RUB) = LUF;
1614  Reflect(RB, RUB) = LUF;
1615  Reflect(RF, RUB) = LUF;
1616 
1617  Reflect(LB, RUF) = LUB;
1618  Reflect(LF, RUF) = LUB;
1619  Reflect(RB, RUF) = LUB;
1620  Reflect(RF, RUF) = LUB;
1621 
1622 
1623  Reflect(DB, LDB) = LUF;
1624  Reflect(DF, LDB) = LUF;
1625  Reflect(UB, LDB) = LUF;
1626  Reflect(UF, LDB) = LUF;
1627 
1628  Reflect(DB, LDF) = LUB;
1629  Reflect(DF, LDF) = LUB;
1630  Reflect(UB, LDF) = LUB;
1631  Reflect(UF, LDF) = LUB;
1632 
1633  Reflect(DB, LUB) = LDF;
1634  Reflect(DF, LUB) = LDF;
1635  Reflect(UB, LUB) = LDF;
1636  Reflect(UF, LUB) = LDF;
1637 
1638  Reflect(DB, LUF) = LDB;
1639  Reflect(DF, LUF) = LDB;
1640  Reflect(UB, LUF) = LDB;
1641  Reflect(UF, LUF) = LDB;
1642 
1643  Reflect(DB, RDB) = RUF;
1644  Reflect(DF, RDB) = RUF;
1645  Reflect(UB, RDB) = RUF;
1646  Reflect(UF, RDB) = RUF;
1647 
1648  Reflect(DB, RDF) = RUB;
1649  Reflect(DF, RDF) = RUB;
1650  Reflect(UB, RDF) = RUB;
1651  Reflect(UF, RDF) = RUB;
1652 
1653  Reflect(DB, RUB) = RDF;
1654  Reflect(DF, RUB) = RDF;
1655  Reflect(UB, RUB) = RDF;
1656  Reflect(UF, RUB) = RDF;
1657 
1658  Reflect(DB, RUF) = RDB;
1659  Reflect(DF, RUF) = RDB;
1660  Reflect(UB, RUF) = RDB;
1661  Reflect(UF, RUF) = RDB;
1662 
1663 
1664  // S_base(i,direction) etc.: Initial value/increment for coordinate s[i] on
1665  // the face indicated by direction (L/R/U/D/F/B)
1666  S_base.resize(3, 27);
1667  S_steplo.resize(3, 27);
1668  S_stephi.resize(3, 27);
1669 
1670  S_base(0, L) = -1.0;
1671  S_base(1, L) = -1.0;
1672  S_base(2, L) = -1.0;
1673  S_steplo(0, L) = 0.0;
1674  S_steplo(1, L) = 2.0;
1675  S_steplo(2, L) = 0.0;
1676  S_stephi(0, L) = 0.0;
1677  S_stephi(1, L) = 0.0;
1678  S_stephi(2, L) = 2.0;
1679 
1680  S_base(0, R) = 1.0;
1681  S_base(1, R) = -1.0;
1682  S_base(2, R) = -1.0;
1683  S_steplo(0, R) = 0.0;
1684  S_steplo(1, R) = 2.0;
1685  S_steplo(2, R) = 0.0;
1686  S_stephi(0, R) = 0.0;
1687  S_stephi(1, R) = 0.0;
1688  S_stephi(2, R) = 2.0;
1689 
1690  S_base(0, U) = -1.0;
1691  S_base(1, U) = 1.0;
1692  S_base(2, U) = -1.0;
1693  S_steplo(0, U) = 2.0;
1694  S_steplo(1, U) = 0.0;
1695  S_steplo(2, U) = 0.0;
1696  S_stephi(0, U) = 0.0;
1697  S_stephi(1, U) = 0.0;
1698  S_stephi(2, U) = 2.0;
1699 
1700  S_base(0, D) = -1.0;
1701  S_base(1, D) = -1.0;
1702  S_base(2, D) = -1.0;
1703  S_steplo(0, D) = 2.0;
1704  S_steplo(1, D) = 0.0;
1705  S_steplo(2, D) = 0.0;
1706  S_stephi(0, D) = 0.0;
1707  S_stephi(1, D) = 0.0;
1708  S_stephi(2, D) = 2.0;
1709 
1710  S_base(0, F) = -1.0;
1711  S_base(1, F) = -1.0;
1712  S_base(2, F) = 1.0;
1713  S_steplo(0, F) = 2.0;
1714  S_steplo(1, F) = 0.0;
1715  S_steplo(2, F) = 0.0;
1716  S_stephi(0, F) = 0.0;
1717  S_stephi(1, F) = 2.0;
1718  S_stephi(2, F) = 0.0;
1719 
1720  S_base(0, B) = -1.0;
1721  S_base(1, B) = -1.0;
1722  S_base(2, B) = -1.0;
1723  S_steplo(0, B) = 2.0;
1724  S_steplo(1, B) = 0.0;
1725  S_steplo(2, B) = 0.0;
1726  S_stephi(0, B) = 0.0;
1727  S_stephi(1, B) = 2.0;
1728  S_stephi(2, B) = 0.0;
1729 
1730 
1731  // Relative to the left/down/back vertex in any (father) octree, the
1732  // corresponding vertex in the son specified by \c son_octant has an offset.
1733  // If we project the son_octant's left/down/back vertex onto the
1734  // father's face \c face, it is located at the in-face coordinate
1735  // \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
1736  // \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
1737  S_directlo.resize(27, 27);
1738  S_directhi.resize(27, 27);
1739 
1740  S_directlo(L, LDB) = 0.0;
1741  S_directlo(R, LDB) = 0.0;
1742  S_directlo(U, LDB) = 0.0;
1743  S_directlo(D, LDB) = 0.0;
1744  S_directlo(F, LDB) = 0.0;
1745  S_directlo(B, LDB) = 0.0;
1746 
1747  S_directlo(L, LDF) = 0.0;
1748  S_directlo(R, LDF) = 0.0;
1749  S_directlo(U, LDF) = 0.0;
1750  S_directlo(D, LDF) = 0.0;
1751  S_directlo(F, LDF) = 0.0;
1752  S_directlo(B, LDF) = 0.0;
1753 
1754  S_directlo(L, LUB) = 1.0;
1755  S_directlo(R, LUB) = 1.0;
1756  S_directlo(U, LUB) = 0.0;
1757  S_directlo(D, LUB) = 0.0;
1758  S_directlo(F, LUB) = 0.0;
1759  S_directlo(B, LUB) = 0.0;
1760 
1761  S_directlo(L, LUF) = 1.0;
1762  S_directlo(R, LUF) = 1.0;
1763  S_directlo(U, LUF) = 0.0;
1764  S_directlo(D, LUF) = 0.0;
1765  S_directlo(F, LUF) = 0.0;
1766  S_directlo(B, LUF) = 0.0;
1767 
1768  S_directlo(L, RDB) = 0.0;
1769  S_directlo(R, RDB) = 0.0;
1770  S_directlo(U, RDB) = 1.0;
1771  S_directlo(D, RDB) = 1.0;
1772  S_directlo(F, RDB) = 1.0;
1773  S_directlo(B, RDB) = 1.0;
1774 
1775  S_directlo(L, RDF) = 0.0;
1776  S_directlo(R, RDF) = 0.0;
1777  S_directlo(U, RDF) = 1.0;
1778  S_directlo(D, RDF) = 1.0;
1779  S_directlo(F, RDF) = 1.0;
1780  S_directlo(B, RDF) = 1.0;
1781 
1782  S_directlo(L, RUB) = 1.0;
1783  S_directlo(R, RUB) = 1.0;
1784  S_directlo(U, RUB) = 1.0;
1785  S_directlo(D, RUB) = 1.0;
1786  S_directlo(F, RUB) = 1.0;
1787  S_directlo(B, RUB) = 1.0;
1788 
1789  S_directlo(L, RUF) = 1.0;
1790  S_directlo(R, RUF) = 1.0;
1791  S_directlo(U, RUF) = 1.0;
1792  S_directlo(D, RUF) = 1.0;
1793  S_directlo(F, RUF) = 1.0;
1794  S_directlo(B, RUF) = 1.0;
1795 
1796 
1797  S_directhi(L, LDB) = 0.0;
1798  S_directhi(R, LDB) = 0.0;
1799  S_directhi(U, LDB) = 0.0;
1800  S_directhi(D, LDB) = 0.0;
1801  S_directhi(F, LDB) = 0.0;
1802  S_directhi(B, LDB) = 0.0;
1803 
1804  S_directhi(L, LDF) = 1.0;
1805  S_directhi(R, LDF) = 1.0;
1806  S_directhi(U, LDF) = 1.0;
1807  S_directhi(D, LDF) = 1.0;
1808  S_directhi(F, LDF) = 0.0;
1809  S_directhi(B, LDF) = 0.0;
1810 
1811  S_directhi(L, LUB) = 0.0;
1812  S_directhi(R, LUB) = 0.0;
1813  S_directhi(U, LUB) = 0.0;
1814  S_directhi(D, LUB) = 0.0;
1815  S_directhi(F, LUB) = 1.0;
1816  S_directhi(B, LUB) = 1.0;
1817 
1818  S_directhi(L, LUF) = 1.0;
1819  S_directhi(R, LUF) = 1.0;
1820  S_directhi(U, LUF) = 1.0;
1821  S_directhi(D, LUF) = 1.0;
1822  S_directhi(F, LUF) = 1.0;
1823  S_directhi(B, LUF) = 1.0;
1824 
1825  S_directhi(L, RDB) = 0.0;
1826  S_directhi(R, RDB) = 0.0;
1827  S_directhi(U, RDB) = 0.0;
1828  S_directhi(D, RDB) = 0.0;
1829  S_directhi(F, RDB) = 0.0;
1830  S_directhi(B, RDB) = 0.0;
1831 
1832  S_directhi(L, RDF) = 1.0;
1833  S_directhi(R, RDF) = 1.0;
1834  S_directhi(U, RDF) = 1.0;
1835  S_directhi(D, RDF) = 1.0;
1836  S_directhi(F, RDF) = 0.0;
1837  S_directhi(B, RDF) = 0.0;
1838 
1839  S_directhi(L, RUB) = 0.0;
1840  S_directhi(R, RUB) = 0.0;
1841  S_directhi(U, RUB) = 0.0;
1842  S_directhi(D, RUB) = 0.0;
1843  S_directhi(F, RUB) = 1.0;
1844  S_directhi(B, RUB) = 1.0;
1845 
1846  S_directhi(L, RUF) = 1.0;
1847  S_directhi(R, RUF) = 1.0;
1848  S_directhi(U, RUF) = 1.0;
1849  S_directhi(D, RUF) = 1.0;
1850  S_directhi(F, RUF) = 1.0;
1851  S_directhi(B, RUF) = 1.0;
1852 
1853 
1854  // S_base_edge(i,direction): Initial value/increment for coordinate s[i] on
1855  // the edge indicated by direction (LB,RB,...)
1856  S_base_edge.resize(3, 27);
1857  S_step_edge.resize(3, 27);
1858 
1859  S_base_edge(0, LB) = -1.0;
1860  S_base_edge(1, LB) = -1.0;
1861  S_base_edge(2, LB) = -1.0;
1862  S_step_edge(0, LB) = 0.0;
1863  S_step_edge(1, LB) = 2.0;
1864  S_step_edge(2, LB) = 0.0;
1865 
1866  S_base_edge(0, RB) = 1.0;
1867  S_base_edge(1, RB) = -1.0;
1868  S_base_edge(2, RB) = -1.0;
1869  S_step_edge(0, RB) = 0.0;
1870  S_step_edge(1, RB) = 2.0;
1871  S_step_edge(2, RB) = 0.0;
1872 
1873  S_base_edge(0, DB) = -1.0;
1874  S_base_edge(1, DB) = -1.0;
1875  S_base_edge(2, DB) = -1.0;
1876  S_step_edge(0, DB) = 2.0;
1877  S_step_edge(1, DB) = 0.0;
1878  S_step_edge(2, DB) = 0.0;
1879 
1880  S_base_edge(0, UB) = -1.0;
1881  S_base_edge(1, UB) = 1.0;
1882  S_base_edge(2, UB) = -1.0;
1883  S_step_edge(0, UB) = 2.0;
1884  S_step_edge(1, UB) = 0.0;
1885  S_step_edge(2, UB) = 0.0;
1886 
1887  S_base_edge(0, LD) = -1.0;
1888  S_base_edge(1, LD) = -1.0;
1889  S_base_edge(2, LD) = -1.0;
1890  S_step_edge(0, LD) = 0.0;
1891  S_step_edge(1, LD) = 0.0;
1892  S_step_edge(2, LD) = 2.0;
1893 
1894  S_base_edge(0, RD) = 1.0;
1895  S_base_edge(1, RD) = -1.0;
1896  S_base_edge(2, RD) = -1.0;
1897  S_step_edge(0, RD) = 0.0;
1898  S_step_edge(1, RD) = 0.0;
1899  S_step_edge(2, RD) = 2.0;
1900 
1901  S_base_edge(0, LU) = -1.0;
1902  S_base_edge(1, LU) = 1.0;
1903  S_base_edge(2, LU) = -1.0;
1904  S_step_edge(0, LU) = 0.0;
1905  S_step_edge(1, LU) = 0.0;
1906  S_step_edge(2, LU) = 2.0;
1907 
1908  S_base_edge(0, RU) = 1.0;
1909  S_base_edge(1, RU) = 1.0;
1910  S_base_edge(2, RU) = -1.0;
1911  S_step_edge(0, RU) = 0.0;
1912  S_step_edge(1, RU) = 0.0;
1913  S_step_edge(2, RU) = 2.0;
1914 
1915 
1916  S_base_edge(0, LF) = -1.0;
1917  S_base_edge(1, LF) = -1.0;
1918  S_base_edge(2, LF) = 1.0;
1919  S_step_edge(0, LF) = 0.0;
1920  S_step_edge(1, LF) = 2.0;
1921  S_step_edge(2, LF) = 0.0;
1922 
1923  S_base_edge(0, RF) = 1.0;
1924  S_base_edge(1, RF) = -1.0;
1925  S_base_edge(2, RF) = 1.0;
1926  S_step_edge(0, RF) = 0.0;
1927  S_step_edge(1, RF) = 2.0;
1928  S_step_edge(2, RF) = 0.0;
1929 
1930  S_base_edge(0, DF) = -1.0;
1931  S_base_edge(1, DF) = -1.0;
1932  S_base_edge(2, DF) = 1.0;
1933  S_step_edge(0, DF) = 2.0;
1934  S_step_edge(1, DF) = 0.0;
1935  S_step_edge(2, DF) = 0.0;
1936 
1937  S_base_edge(0, UF) = -1.0;
1938  S_base_edge(1, UF) = 1.0;
1939  S_base_edge(2, UF) = 1.0;
1940  S_step_edge(0, UF) = 2.0;
1941  S_step_edge(1, UF) = 0.0;
1942  S_step_edge(2, UF) = 0.0;
1943 
1944 
1945  // Relative to the left/down/back vertex in any (father) octree, the
1946  // corresponding vertex in the son specified by \c son_octant has an offset.
1947  // If we project the son_octant's left/down/back vertex onto the
1948  // father's edge \c edge, it is located at the in-face coordinate
1949  // \c s_lo = h/2 \c S_direct_edge(edge,son_octant).
1950  S_direct_edge.resize(27, 27);
1951  S_direct_edge(LB, LDB) = 0.0;
1952  S_direct_edge(RB, LDB) = 0.0;
1953  S_direct_edge(DB, LDB) = 0.0;
1954  S_direct_edge(UB, LDB) = 0.0;
1955  S_direct_edge(LD, LDB) = 0.0;
1956  S_direct_edge(RD, LDB) = 0.0;
1957  S_direct_edge(LU, LDB) = 0.0;
1958  S_direct_edge(RU, LDB) = 0.0;
1959  S_direct_edge(LF, LDB) = 0.0;
1960  S_direct_edge(RF, LDB) = 0.0;
1961  S_direct_edge(DF, LDB) = 0.0;
1962  S_direct_edge(UF, LDB) = 0.0;
1963 
1964  S_direct_edge(LB, RDB) = 0.0;
1965  S_direct_edge(RB, RDB) = 0.0;
1966  S_direct_edge(DB, RDB) = 1.0;
1967  S_direct_edge(UB, RDB) = 1.0;
1968  S_direct_edge(LD, RDB) = 0.0;
1969  S_direct_edge(RD, RDB) = 0.0;
1970  S_direct_edge(LU, RDB) = 0.0;
1971  S_direct_edge(RU, RDB) = 0.0;
1972  S_direct_edge(LF, RDB) = 0.0;
1973  S_direct_edge(RF, RDB) = 0.0;
1974  S_direct_edge(DF, RDB) = 1.0;
1975  S_direct_edge(UF, RDB) = 1.0;
1976 
1977  S_direct_edge(LB, LUB) = 1.0;
1978  S_direct_edge(RB, LUB) = 1.0;
1979  S_direct_edge(DB, LUB) = 0.0;
1980  S_direct_edge(UB, LUB) = 0.0;
1981  S_direct_edge(LD, LUB) = 0.0;
1982  S_direct_edge(RD, LUB) = 0.0;
1983  S_direct_edge(LU, LUB) = 0.0;
1984  S_direct_edge(RU, LUB) = 0.0;
1985  S_direct_edge(LF, LUB) = 1.0;
1986  S_direct_edge(RF, LUB) = 1.0;
1987  S_direct_edge(DF, LUB) = 0.0;
1988  S_direct_edge(UF, LUB) = 0.0;
1989 
1990  S_direct_edge(LB, RUB) = 1.0;
1991  S_direct_edge(RB, RUB) = 1.0;
1992  S_direct_edge(DB, RUB) = 1.0;
1993  S_direct_edge(UB, RUB) = 1.0;
1994  S_direct_edge(LD, RUB) = 0.0;
1995  S_direct_edge(RD, RUB) = 0.0;
1996  S_direct_edge(LU, RUB) = 0.0;
1997  S_direct_edge(RU, RUB) = 0.0;
1998  S_direct_edge(LF, RUB) = 1.0;
1999  S_direct_edge(RF, RUB) = 1.0;
2000  S_direct_edge(DF, RUB) = 1.0;
2001  S_direct_edge(UF, RUB) = 1.0;
2002 
2003 
2004  S_direct_edge(LB, LDF) = 0.0;
2005  S_direct_edge(RB, LDF) = 0.0;
2006  S_direct_edge(DB, LDF) = 0.0;
2007  S_direct_edge(UB, LDF) = 0.0;
2008  S_direct_edge(LD, LDF) = 1.0;
2009  S_direct_edge(RD, LDF) = 1.0;
2010  S_direct_edge(LU, LDF) = 1.0;
2011  S_direct_edge(RU, LDF) = 1.0;
2012  S_direct_edge(LF, LDF) = 0.0;
2013  S_direct_edge(RF, LDF) = 0.0;
2014  S_direct_edge(DF, LDF) = 0.0;
2015  S_direct_edge(UF, LDF) = 0.0;
2016 
2017  S_direct_edge(LB, RDF) = 0.0;
2018  S_direct_edge(RB, RDF) = 0.0;
2019  S_direct_edge(DB, RDF) = 1.0;
2020  S_direct_edge(UB, RDF) = 1.0;
2021  S_direct_edge(LD, RDF) = 1.0;
2022  S_direct_edge(RD, RDF) = 1.0;
2023  S_direct_edge(LU, RDF) = 1.0;
2024  S_direct_edge(RU, RDF) = 1.0;
2025  S_direct_edge(LF, RDF) = 0.0;
2026  S_direct_edge(RF, RDF) = 0.0;
2027  S_direct_edge(DF, RDF) = 1.0;
2028  S_direct_edge(UF, RDF) = 1.0;
2029 
2030  S_direct_edge(LB, LUF) = 1.0;
2031  S_direct_edge(RB, LUF) = 1.0;
2032  S_direct_edge(DB, LUF) = 0.0;
2033  S_direct_edge(UB, LUF) = 0.0;
2034  S_direct_edge(LD, LUF) = 1.0;
2035  S_direct_edge(RD, LUF) = 1.0;
2036  S_direct_edge(LU, LUF) = 1.0;
2037  S_direct_edge(RU, LUF) = 1.0;
2038  S_direct_edge(LF, LUF) = 1.0;
2039  S_direct_edge(RF, LUF) = 1.0;
2040  S_direct_edge(DF, LUF) = 0.0;
2041  S_direct_edge(UF, LUF) = 0.0;
2042 
2043  S_direct_edge(LB, RUF) = 1.0;
2044  S_direct_edge(RB, RUF) = 1.0;
2045  S_direct_edge(DB, RUF) = 1.0;
2046  S_direct_edge(UB, RUF) = 1.0;
2047  S_direct_edge(LD, RUF) = 1.0;
2048  S_direct_edge(RD, RUF) = 1.0;
2049  S_direct_edge(LU, RUF) = 1.0;
2050  S_direct_edge(RU, RUF) = 1.0;
2051  S_direct_edge(LF, RUF) = 1.0;
2052  S_direct_edge(RF, RUF) = 1.0;
2053  S_direct_edge(DF, RUF) = 1.0;
2054  S_direct_edge(UF, RUF) = 1.0;
2055 
2056 
2057  // Public static data:
2058  //-------------------
2059 
2060  // Translate (enumerated) directions into strings
2061  Direct_string.resize(27);
2062  Direct_string[LDB] = "LDB";
2063  Direct_string[LDF] = "LDF";
2064  Direct_string[LUB] = "LUB";
2065  Direct_string[LUF] = "LUF";
2066  Direct_string[RDB] = "RDB";
2067  Direct_string[RDF] = "RDF";
2068  Direct_string[RUB] = "RUB";
2069  Direct_string[RUF] = "RUF";
2070 
2071 
2072  Direct_string[L] = "L";
2073  Direct_string[R] = "R";
2074  Direct_string[U] = "U";
2075  Direct_string[D] = "D";
2076  Direct_string[F] = "F";
2077  Direct_string[B] = "B";
2078 
2079  Direct_string[LU] = "LU";
2080  Direct_string[LD] = "LD";
2081  Direct_string[LF] = "LF";
2082  Direct_string[LB] = "LB";
2083  Direct_string[RU] = "RU";
2084  Direct_string[RD] = "RD";
2085  Direct_string[RF] = "RF";
2086  Direct_string[RB] = "RB";
2087  Direct_string[UF] = "UF";
2088  Direct_string[UB] = "UB";
2089  Direct_string[DF] = "DF";
2090  Direct_string[DB] = "DB";
2091 
2092  Direct_string[OMEGA] = "OMEGA";
2093 
2094 
2095  // Get opposite face, e.g. Reflect_face(L)=R
2096  Reflect_face.resize(27);
2097  Reflect_face[L] = R;
2098  Reflect_face[R] = L;
2099  Reflect_face[U] = D;
2100  Reflect_face[D] = U;
2101  Reflect_face[B] = F;
2102  Reflect_face[F] = B;
2103 
2104  // Get opposite edge, e.g. Reflect_edge(DB)=UF
2105  Reflect_edge.resize(27);
2106  Reflect_edge[LB] = RF;
2107  Reflect_edge[RB] = LF;
2108  Reflect_edge[DB] = UF;
2109  Reflect_edge[UB] = DF;
2110 
2111  Reflect_edge[LD] = RU;
2112  Reflect_edge[RD] = LU;
2113  Reflect_edge[LU] = RD;
2114  Reflect_edge[RU] = LD;
2115 
2116  Reflect_edge[LF] = RB;
2117  Reflect_edge[RF] = LB;
2118  Reflect_edge[DF] = UB;
2119  Reflect_edge[UF] = DB;
2120 
2121  // Get opposite vertex, e.g. Reflect_vertex(LDB)=RUF
2122  Reflect_vertex.resize(27);
2123  Reflect_vertex[LDB] = RUF;
2124  Reflect_vertex[RUF] = LDB;
2125  Reflect_vertex[RDB] = LUF;
2126  Reflect_vertex[LUF] = RDB;
2127  Reflect_vertex[LUB] = RDF;
2128  Reflect_vertex[RDF] = LUB;
2129  Reflect_vertex[RUB] = LDF;
2130  Reflect_vertex[LDF] = RUB;
2131 
2132 
2133  // Vertices at ends of edges
2134  Vertex_at_end_of_edge.resize(27);
2135 
2136  Vertex_at_end_of_edge[DB].resize(2);
2137  Vertex_at_end_of_edge[DB][0] = LDB; // Pattern: both other indices
2138  Vertex_at_end_of_edge[DB][1] = RDB;
2139 
2140  Vertex_at_end_of_edge[UB].resize(2);
2141  Vertex_at_end_of_edge[UB][0] = LUB; // Pattern: both other indices
2142  Vertex_at_end_of_edge[UB][1] = RUB;
2143 
2144 
2145  Vertex_at_end_of_edge[LB].resize(2);
2146  Vertex_at_end_of_edge[LB][0] = LUB; // Pattern: both other indices
2147  Vertex_at_end_of_edge[LB][1] = LDB;
2148 
2149  Vertex_at_end_of_edge[RB].resize(2);
2150  Vertex_at_end_of_edge[RB][0] = RUB; // Pattern: both other indices
2151  Vertex_at_end_of_edge[RB][1] = RDB;
2152 
2153 
2154  Vertex_at_end_of_edge[LD].resize(2);
2155  Vertex_at_end_of_edge[LD][0] = LDF; // Pattern: both other indices
2156  Vertex_at_end_of_edge[LD][1] = LDB;
2157 
2158  Vertex_at_end_of_edge[RD].resize(2);
2159  Vertex_at_end_of_edge[RD][0] = RDF; // Pattern: both other indices
2160  Vertex_at_end_of_edge[RD][1] = RDB;
2161 
2162  Vertex_at_end_of_edge[LU].resize(2);
2163  Vertex_at_end_of_edge[LU][0] = LUF; // Pattern: both other indices
2164  Vertex_at_end_of_edge[LU][1] = LUB;
2165 
2166  Vertex_at_end_of_edge[RU].resize(2);
2167  Vertex_at_end_of_edge[RU][0] = RUF; // Pattern: both other indices
2168  Vertex_at_end_of_edge[RU][1] = RUB;
2169 
2170 
2171  Vertex_at_end_of_edge[DF].resize(2);
2172  Vertex_at_end_of_edge[DF][0] = LDF; // Pattern: both other indices
2173  Vertex_at_end_of_edge[DF][1] = RDF;
2174 
2175  Vertex_at_end_of_edge[UF].resize(2);
2176  Vertex_at_end_of_edge[UF][0] = LUF; // Pattern: both other indices
2177  Vertex_at_end_of_edge[UF][1] = RUF;
2178 
2179  Vertex_at_end_of_edge[LF].resize(2);
2180  Vertex_at_end_of_edge[LF][0] = LUF; // Pattern: both other indices
2181  Vertex_at_end_of_edge[LF][1] = LDF;
2182 
2183  Vertex_at_end_of_edge[RF].resize(2);
2184  Vertex_at_end_of_edge[RF][0] = RUF; // Pattern: both other indices
2185  Vertex_at_end_of_edge[RF][1] = RDF;
2186 
2187 
2188  // Initialisation of the values of Vector_to_direction
2189  Vector<int> vect(3);
2190  int elem;
2191 
2192  for (int i = -1; i < 2; i++)
2193  {
2194  for (int j = -1; j < 2; j++)
2195  {
2196  for (int k = -1; k < 2; k++)
2197  {
2198  vect[0] = i;
2199  vect[1] = j;
2200  vect[2] = k;
2201  int num_elem = 0;
2202 
2203  // To put a number on the vector (i,j,k), we assume that that
2204  // the vector (i+1,j+1,k+1) represents the decomposition
2205  // of the number of the corresponding direction in base 3.
2206  num_elem = (i + 1) * 9 + (j + 1) * 3 + (k + 1);
2207 
2208  // for each number we have the corresponding element
2209  switch (num_elem)
2210  {
2211  case 6:
2212  elem = LUB;
2213  break;
2214  case 24:
2215  elem = RUB;
2216  break;
2217  case 26:
2218  elem = RUF;
2219  break;
2220  case 8:
2221  elem = LUF;
2222  break;
2223  case 0:
2224  elem = LDB;
2225  break;
2226  case 18:
2227  elem = RDB;
2228  break;
2229  case 20:
2230  elem = RDF;
2231  break;
2232  case 2:
2233  elem = LDF;
2234  break;
2235  case 25:
2236  elem = RU;
2237  break;
2238  case 23:
2239  elem = RF;
2240  break;
2241  case 19:
2242  elem = RD;
2243  break;
2244  case 21:
2245  elem = RB;
2246  break;
2247  case 7:
2248  elem = LU;
2249  break;
2250  case 5:
2251  elem = LF;
2252  break;
2253  case 1:
2254  elem = LD;
2255  break;
2256  case 3:
2257  elem = LB;
2258  break;
2259  case 17:
2260  elem = UF;
2261  break;
2262  case 15:
2263  elem = UB;
2264  break;
2265  case 11:
2266  elem = DF;
2267  break;
2268  case 9:
2269  elem = DB;
2270  break;
2271  case 16:
2272  elem = U;
2273  break;
2274  case 10:
2275  elem = D;
2276  break;
2277  case 22:
2278  elem = R;
2279  break;
2280  case 4:
2281  elem = L;
2282  break;
2283  case 14:
2284  elem = F;
2285  break;
2286  case 12:
2287  elem = B;
2288  break;
2289  case 13:
2290  elem = OMEGA;
2291  break;
2292  default:
2293  elem = OMEGA;
2294  oomph_info << "there might be a problem with Vector_to_direction"
2295  << std::endl;
2296  break;
2297  }
2298  Vector_to_direction[vect] = elem;
2299  }
2300  }
2301  }
2302 
2303 
2304  // Initialisation of Direction_to_vector:
2305  // Translate Octant, face, edge into direction vector using the
2306  // value of Direct_string; Direction_to_vector[U]={0,1,0};
2307  Direction_to_vector.resize(27);
2308  for (int i = LDB; i <= F; i++)
2309  {
2310  Direction_to_vector[i].resize(3);
2311  // Initialisation to 0;
2312  for (int j = 0; j < 3; j++)
2313  {
2314  Direction_to_vector[i][j] = 0;
2315  }
2316 
2317  // Use +1 or -1 to indicate the relevant components of
2318  // the vector Direction
2319  for (unsigned j = 0; j < 3; j++)
2320  {
2321  if (Direct_string[i].length() > j)
2322  {
2323  switch (Direct_string[i][j])
2324  {
2325  case 'R':
2326  Direction_to_vector[i][0] = 1;
2327  break;
2328  case 'L':
2329  Direction_to_vector[i][0] = -1;
2330  break;
2331  case 'U':
2332  Direction_to_vector[i][1] = 1;
2333  break;
2334  case 'D':
2335  Direction_to_vector[i][1] = -1;
2336  break;
2337  case 'F':
2338  Direction_to_vector[i][2] = 1;
2339  break;
2340  case 'B':
2341  Direction_to_vector[i][2] = -1;
2342  break;
2343  default:
2344  oomph_info << "Direction Error !!" << std::endl;
2345  }
2346  }
2347  }
2348  }
2349 
2350 
2351  // Setup map that works out required rotations based on
2352  //-----------------------------------------------------
2353  // adjacent edge vertices
2354  //-----------------------
2355 
2356  int new_up, new_right;
2357  int new_vertex;
2358 
2359 
2360  // Map that stores the set of rotations (as a pairs consisting of
2361  // the up_equivalent and the right_equivalent) that move the vertex
2362  // specified by the first entry in key pair to the position of the second
2363  // one:
2364  std::map<std::pair<int, int>, std::set<std::pair<int, int>>>
2365  required_rotation;
2366 
2367  // Loop over all vertices
2368  for (int vertex = LDB; vertex <= RUF; vertex++)
2369  {
2370  new_up = U;
2371  new_right = R;
2372  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2373  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2374  std::make_pair(new_up, new_right));
2375 
2376  new_up = U;
2377  new_right = B;
2378  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2379  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2380  std::make_pair(new_up, new_right));
2381 
2382  new_up = U;
2383  new_right = L;
2384  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2385  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2386  std::make_pair(new_up, new_right));
2387 
2388  new_up = U;
2389  new_right = F;
2390  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2391  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2392  std::make_pair(new_up, new_right));
2393 
2394 
2395  new_up = D;
2396  new_right = R;
2397  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2398  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2399  std::make_pair(new_up, new_right));
2400 
2401  new_up = D;
2402  new_right = B;
2403  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2404  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2405  std::make_pair(new_up, new_right));
2406 
2407  new_up = D;
2408  new_right = L;
2409  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2410  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2411  std::make_pair(new_up, new_right));
2412 
2413  new_up = D;
2414  new_right = F;
2415  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2416  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2417  std::make_pair(new_up, new_right));
2418 
2419 
2420  new_up = R;
2421  new_right = D;
2422  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2423  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2424  std::make_pair(new_up, new_right));
2425 
2426  new_up = R;
2427  new_right = B;
2428  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2429  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2430  std::make_pair(new_up, new_right));
2431 
2432  new_up = R;
2433  new_right = U;
2434  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2435  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2436  std::make_pair(new_up, new_right));
2437 
2438  new_up = R;
2439  new_right = F;
2440  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2441  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2442  std::make_pair(new_up, new_right));
2443 
2444 
2445  new_up = L;
2446  new_right = D;
2447  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2448  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2449  std::make_pair(new_up, new_right));
2450 
2451  new_up = L;
2452  new_right = B;
2453  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2454  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2455  std::make_pair(new_up, new_right));
2456 
2457  new_up = L;
2458  new_right = U;
2459  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2460  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2461  std::make_pair(new_up, new_right));
2462 
2463  new_up = L;
2464  new_right = F;
2465  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2466  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2467  std::make_pair(new_up, new_right));
2468 
2469 
2470  new_up = F;
2471  new_right = R;
2472  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2473  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2474  std::make_pair(new_up, new_right));
2475 
2476  new_up = F;
2477  new_right = L;
2478  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2479  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2480  std::make_pair(new_up, new_right));
2481 
2482  new_up = F;
2483  new_right = U;
2484  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2485  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2486  std::make_pair(new_up, new_right));
2487 
2488  new_up = F;
2489  new_right = D;
2490  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2491  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2492  std::make_pair(new_up, new_right));
2493 
2494 
2495  new_up = B;
2496  new_right = R;
2497  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2498  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2499  std::make_pair(new_up, new_right));
2500 
2501  new_up = B;
2502  new_right = L;
2503  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2504  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2505  std::make_pair(new_up, new_right));
2506 
2507  new_up = B;
2508  new_right = U;
2509  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2510  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2511  std::make_pair(new_up, new_right));
2512 
2513  new_up = B;
2514  new_right = D;
2515  new_vertex = OcTree::rotate(new_up, new_right, vertex);
2516  required_rotation[std::make_pair(vertex, new_vertex)].insert(
2517  std::make_pair(new_up, new_right));
2518  }
2519 
2520 
2521  // Each vertex is part of three edges. This container stores the
2522  // vertices in each of the three edge neighbours that are
2523  // adjacent to this node if there's no relative rotation between
2524  // the elements.
2525  std::map<int, Vector<int>> vertex_in_edge_neighbour;
2526 
2527 
2528  // Each vertex is part of three edges. This container stores the
2529  // vertices at the other end of the edge
2530  std::map<int, Vector<int>> other_vertex_on_edge;
2531 
2532  // Each vertex is part of three edges. This container stores the
2533  // vertices in the adjacent element at the other end of the edge
2534  // assuming there are no rotations between the elements.
2535  std::map<int, Vector<int>> other_vertex_in_edge_neighbour;
2536 
2537 
2538  vertex_in_edge_neighbour[LDB].resize(3);
2539  vertex_in_edge_neighbour[LDB][0] =
2540  LUF; // Pattern: exactly one letter matches
2541  vertex_in_edge_neighbour[LDB][1] = RDF;
2542  vertex_in_edge_neighbour[LDB][2] = RUB;
2543 
2544  other_vertex_on_edge[LDB].resize(3);
2545  other_vertex_on_edge[LDB][0] =
2546  RDB; // Pattern: opposite of the matching letter
2547  other_vertex_on_edge[LDB][1] = LUB;
2548  other_vertex_on_edge[LDB][2] = LDF;
2549 
2550  other_vertex_in_edge_neighbour[LDB].resize(3);
2551  other_vertex_in_edge_neighbour[LDB][0] = RUF; // Pattern: full reflection
2552  other_vertex_in_edge_neighbour[LDB][1] = RUF;
2553  other_vertex_in_edge_neighbour[LDB][2] = RUF;
2554 
2555 
2556  vertex_in_edge_neighbour[RDB].resize(3);
2557  vertex_in_edge_neighbour[RDB][0] =
2558  RUF; // Pattern: exactly one letter matches
2559  vertex_in_edge_neighbour[RDB][1] = LDF;
2560  vertex_in_edge_neighbour[RDB][2] = LUB;
2561 
2562  other_vertex_on_edge[RDB].resize(3);
2563  other_vertex_on_edge[RDB][0] =
2564  LDB; // Pattern: opposite of the matching letter
2565  other_vertex_on_edge[RDB][1] = RUB;
2566  other_vertex_on_edge[RDB][2] = RDF;
2567 
2568  other_vertex_in_edge_neighbour[RDB].resize(3);
2569  other_vertex_in_edge_neighbour[RDB][0] = LUF; // Pattern: full reflection
2570  other_vertex_in_edge_neighbour[RDB][1] = LUF;
2571  other_vertex_in_edge_neighbour[RDB][2] = LUF;
2572 
2573 
2574  vertex_in_edge_neighbour[LUB].resize(3);
2575  vertex_in_edge_neighbour[LUB][0] =
2576  LDF; // Pattern: exactly one letter matches
2577  vertex_in_edge_neighbour[LUB][1] = RUF;
2578  vertex_in_edge_neighbour[LUB][2] = RDB;
2579 
2580  other_vertex_on_edge[LUB].resize(3);
2581  other_vertex_on_edge[LUB][0] =
2582  RUB; // Pattern: opposite of the matching letter
2583  other_vertex_on_edge[LUB][1] = LDB;
2584  other_vertex_on_edge[LUB][2] = LUF;
2585 
2586  other_vertex_in_edge_neighbour[LUB].resize(3);
2587  other_vertex_in_edge_neighbour[LUB][0] = RDF; // Pattern: full reflection
2588  other_vertex_in_edge_neighbour[LUB][1] = RDF;
2589  other_vertex_in_edge_neighbour[LUB][2] = RDF;
2590 
2591 
2592  vertex_in_edge_neighbour[RUB].resize(3);
2593  vertex_in_edge_neighbour[RUB][0] =
2594  RDF; // Pattern: exactly one letter matches
2595  vertex_in_edge_neighbour[RUB][1] = LUF;
2596  vertex_in_edge_neighbour[RUB][2] = LDB;
2597 
2598  other_vertex_on_edge[RUB].resize(3);
2599  other_vertex_on_edge[RUB][0] =
2600  LUB; // Pattern: opposite of the matching letter
2601  other_vertex_on_edge[RUB][1] = RDB;
2602  other_vertex_on_edge[RUB][2] = RUF;
2603 
2604  other_vertex_in_edge_neighbour[RUB].resize(3);
2605  other_vertex_in_edge_neighbour[RUB][0] = LDF; // Pattern: full reflection
2606  other_vertex_in_edge_neighbour[RUB][1] = LDF;
2607  other_vertex_in_edge_neighbour[RUB][2] = LDF;
2608 
2609 
2610  vertex_in_edge_neighbour[LDF].resize(3);
2611  vertex_in_edge_neighbour[LDF][0] =
2612  LUB; // Pattern: exactly one letter matches
2613  vertex_in_edge_neighbour[LDF][1] = RDB;
2614  vertex_in_edge_neighbour[LDF][2] = RUF;
2615 
2616  other_vertex_on_edge[LDF].resize(3);
2617  other_vertex_on_edge[LDF][0] =
2618  RDF; // Pattern: opposite of the matching letter
2619  other_vertex_on_edge[LDF][1] = LUF;
2620  other_vertex_on_edge[LDF][2] = LDB;
2621 
2622  other_vertex_in_edge_neighbour[LDF].resize(3);
2623  other_vertex_in_edge_neighbour[LDF][0] = RUB; // Pattern: full reflection
2624  other_vertex_in_edge_neighbour[LDF][1] = RUB;
2625  other_vertex_in_edge_neighbour[LDF][2] = RUB;
2626 
2627 
2628  vertex_in_edge_neighbour[RDF].resize(3);
2629  vertex_in_edge_neighbour[RDF][0] =
2630  RUB; // Pattern: exactly one letter matches
2631  vertex_in_edge_neighbour[RDF][1] = LDB;
2632  vertex_in_edge_neighbour[RDF][2] = LUF;
2633 
2634  other_vertex_on_edge[RDF].resize(3);
2635  other_vertex_on_edge[RDF][0] =
2636  LDF; // Pattern: opposite of the matching letter
2637  other_vertex_on_edge[RDF][1] = RUF;
2638  other_vertex_on_edge[RDF][2] = RDB;
2639 
2640  other_vertex_in_edge_neighbour[RDF].resize(3);
2641  other_vertex_in_edge_neighbour[RDF][0] = LUB; // Pattern: full reflection
2642  other_vertex_in_edge_neighbour[RDF][1] = LUB;
2643  other_vertex_in_edge_neighbour[RDF][2] = LUB;
2644 
2645 
2646  vertex_in_edge_neighbour[LUF].resize(3);
2647  vertex_in_edge_neighbour[LUF][0] =
2648  LDB; // Pattern: exactly one letter matches
2649  vertex_in_edge_neighbour[LUF][1] = RUB;
2650  vertex_in_edge_neighbour[LUF][2] = RDF;
2651 
2652  other_vertex_on_edge[LUF].resize(3);
2653  other_vertex_on_edge[LUF][0] =
2654  RUF; // Pattern: opposite of the matching letter
2655  other_vertex_on_edge[LUF][1] = LDF;
2656  other_vertex_on_edge[LUF][2] = LUB;
2657 
2658  other_vertex_in_edge_neighbour[LUF].resize(3);
2659  other_vertex_in_edge_neighbour[LUF][0] = RDB; // Pattern: full reflection
2660  other_vertex_in_edge_neighbour[LUF][1] = RDB;
2661  other_vertex_in_edge_neighbour[LUF][2] = RDB;
2662 
2663 
2664  vertex_in_edge_neighbour[RUF].resize(3);
2665  vertex_in_edge_neighbour[RUF][0] =
2666  RDB; // Pattern: exactly one letter matches
2667  vertex_in_edge_neighbour[RUF][1] = LUB;
2668  vertex_in_edge_neighbour[RUF][2] = LDF;
2669 
2670  other_vertex_on_edge[RUF].resize(3);
2671  other_vertex_on_edge[RUF][0] =
2672  LUF; // Pattern: opposite of the matching letter
2673  other_vertex_on_edge[RUF][1] = RDF;
2674  other_vertex_on_edge[RUF][2] = RUB;
2675 
2676  other_vertex_in_edge_neighbour[RUF].resize(3);
2677  other_vertex_in_edge_neighbour[RUF][0] = LDB; // Pattern: full reflection
2678  other_vertex_in_edge_neighbour[RUF][1] = LDB;
2679  other_vertex_in_edge_neighbour[RUF][2] = LDB;
2680 
2681 
2682  // Loop over all vertices in the reference element
2683  for (int vertex = LDB; vertex <= RUF; vertex++)
2684  {
2685  // Loop over the three edges that are connected to this vertex
2686  for (unsigned i = 0; i < 3; i++)
2687  {
2688  // This is the other vertex along this edge
2689  int other_vertex = other_vertex_on_edge[vertex][i];
2690 
2691  // This is the vertex in the edge neighbour that corresponds
2692  // to the present vertex in the reference element (in
2693  // the absence of rotations)
2694  int unrotated_neigh_vertex = vertex_in_edge_neighbour[vertex][i];
2695 
2696  // This is the vertex in the edge neighbour that corresponds
2697  // to the other vertex in the reference element (in
2698  // the absence of rotations)
2699  int unrotated_neigh_other_vertex =
2700  other_vertex_in_edge_neighbour[vertex][i];
2701 
2702  // Loop over all vertices in the neighbour element
2703  for (int neigh_vertex = LDB; neigh_vertex <= RUF; neigh_vertex++)
2704  {
2705  // What rotations would turn the neigh_vertex
2706  // into the unrotated_neigh_vertex?
2707  std::set<std::pair<int, int>> vertex_rot =
2708  required_rotation[std::make_pair(neigh_vertex,
2709  unrotated_neigh_vertex)];
2710 
2711 
2712  // Loop over all "other" vertices in the neighbour element
2713  for (int neigh_other_vertex = LDB; neigh_other_vertex <= RUF;
2714  neigh_other_vertex++)
2715  {
2716  // What rotations would turn the other_neigh_vertex
2717  // into the unrotated_other_neigh_vertex?
2718  std::set<std::pair<int, int>> other_vertex_rot =
2719  required_rotation[std::make_pair(neigh_other_vertex,
2720  unrotated_neigh_other_vertex)];
2721 
2722  // What are the common rotations?
2723  std::set<std::pair<int, int>> common_rotations;
2724 
2725  // Get the intersection of the two sets
2726  std::set_intersection(
2727  vertex_rot.begin(),
2728  vertex_rot.end(),
2729  other_vertex_rot.begin(),
2730  other_vertex_rot.end(),
2731  inserter(common_rotations, common_rotations.begin()));
2732 
2733 
2734  if (common_rotations.size() > 0)
2735  {
2736  for (std::set<std::pair<int, int>>::iterator it =
2737  common_rotations.begin();
2738  it != common_rotations.end();
2739  it++)
2740  {
2741  // Copy into container
2742 
2743  // First: up equivalent:
2744  Up_and_right_equivalent_for_pairs_of_vertices
2745  [std::make_pair(
2746  std::make_pair(vertex, neigh_vertex),
2747  std::make_pair(other_vertex, neigh_other_vertex))]
2748  .first = it->first;
2749 
2750  // Second: Right equivalent
2751  Up_and_right_equivalent_for_pairs_of_vertices
2752  [std::make_pair(
2753  std::make_pair(vertex, neigh_vertex),
2754  std::make_pair(other_vertex, neigh_other_vertex))]
2755  .second = it->second;
2756  }
2757  }
2758  }
2759  }
2760  }
2761  }
2762  }
num_elem
static int num_elem(char *strv, unsigned elem_len, int term_char, int num_term)
Definition: cfortran.h:639
oomph::DenseMatrix::resize
void resize(const unsigned long &n)
Definition: matrices.h:498
oomph::OcTree::Up_and_right_equivalent_for_pairs_of_vertices
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
Definition: octree.h:377
oomph::OcTree::Reflect_vertex
static Vector< int > Reflect_vertex
Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF.
Definition: octree.h:338
oomph::OcTree::Vertex_at_end_of_edge
static Vector< Vector< int > > Vertex_at_end_of_edge
Map of vectors containing the two vertices for each edge.
Definition: octree.h:343
oomph::OcTree::Static_data_has_been_setup
static bool Static_data_has_been_setup
Bool indicating that static member data has been setup.
Definition: octree.h:412
set
void set(Container &c, Position position, const Value &value)
Definition: stdlist_overload.cpp:36

References Colour, Common_face, Cosi, D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, Direct_string, Direction_to_vector, oomph::OcTreeNames::F, i, Is_adjacent, j, k, L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, num_elem(), oomph::Tree::OMEGA, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, Reflect, Reflect_edge, Reflect_face, Reflect_vertex, oomph::DenseMatrix< T >::resize(), oomph::OcTreeNames::RF, rotate(), oomph::OcTreeNames::RU, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, S_base, S_base_edge, S_direct_edge, S_directhi, S_directlo, S_step_edge, S_stephi, S_steplo, set(), Sini, Static_data_has_been_setup, RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, oomph::OcTreeNames::UF, Up_and_right_equivalent_for_pairs_of_vertices, Vector_to_direction, and Vertex_at_end_of_edge.

Referenced by oomph::RefineableBrickMesh< ELEMENT >::RefineableBrickMesh(), and oomph::RefineableEighthSphereMesh< ELEMENT >::RefineableEighthSphereMesh().

◆ vertex_node_to_vector()

Vector< int > oomph::OcTree::vertex_node_to_vector ( const unsigned &  n,
const unsigned &  nnode1d 
)
staticprivate

Returns the vector of the coordinate directions of vertex node number n in an element with nnode1d element per dimension.

This function is used to translate the position of a vertex node (given by his local number n into a vector giving the position of this node in the local coordinates system. It also needs the value of nnode1d to work.

234  {
235 #ifdef PARANOID
236  if ((n != 0) && (n != nnode1d - 1) && (n != (nnode1d - 1) * nnode1d) &&
237  (n != nnode1d * nnode1d - 1) &&
238  (n != (nnode1d * nnode1d) * (nnode1d - 1) + 0) &&
239  (n != (nnode1d * nnode1d) * (nnode1d - 1) + nnode1d - 1) &&
240  (n != (nnode1d * nnode1d) * (nnode1d - 1) + (nnode1d - 1) * nnode1d) &&
241  (n != (nnode1d * nnode1d) * (nnode1d - 1) + nnode1d * nnode1d - 1))
242  {
243  std::ostringstream error_stream;
244  error_stream << "Node " << n
245  << " is not a vertex node in a brick element with "
246  << nnode1d << " nodes along each edge!";
247 
248  throw OomphLibError(
249  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
250  }
251 #endif
252  int a, b, c;
253  Vector<int> result_vect(3);
254  a = n / (nnode1d * nnode1d);
255  b = (n - a * nnode1d * nnode1d) / nnode1d;
256  c = n - a * nnode1d * nnode1d - b * nnode1d;
257 
258  result_vect[0] = 2 * c / (nnode1d - 1) - 1;
259  result_vect[1] = 2 * b / (nnode1d - 1) - 1;
260  result_vect[2] = 2 * a / (nnode1d - 1) - 1;
261 
262  return result_vect;
263  }

References a, b, calibrate::c, n, OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by get_the_other_face().

◆ vertex_to_node_number()

unsigned oomph::OcTree::vertex_to_node_number ( const int &  vertex,
const unsigned &  nnode1d 
)
static

Return the local node number of given vertex [LDB,RDB,...] in an element with nnode1d nodes in each coordinate direction

Return the local node number of given vertex in an element with nnode1d nodes in each coordinate direction.

416  {
417  using namespace OcTreeNames;
418 
419 #ifdef PARANOID
420  if ((vertex != LDB) && (vertex != RDB) && (vertex != LUB) &&
421  (vertex != RUB) && (vertex != LDF) && (vertex != RDF) &&
422  (vertex != LUF) && (vertex != RUF))
423  {
424  std::ostringstream error_stream;
425  error_stream << "Wrong vertex: " << Direct_string[vertex] << std::endl;
426  throw OomphLibError(
427  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
428  }
429 #endif
430 
431  switch (vertex)
432  {
433  case LDB:
434  return 0;
435  break;
436 
437  case RDB:
438  return nnode1d - 1;
439  break;
440 
441 
442  case LUB:
443  return nnode1d * (nnode1d - 1);
444  break;
445 
446  case RUB:
447  return nnode1d * nnode1d - 1;
448  break;
449 
450 
451  case LDF:
452  return nnode1d * nnode1d * (nnode1d - 1);
453  break;
454 
455 
456  case RDF:
457  return (nnode1d * nnode1d + 1) * (nnode1d - 1);
458  break;
459 
460  case LUF:
461  return nnode1d * nnode1d * nnode1d - nnode1d;
462  break;
463 
464  case RUF:
465  return nnode1d * nnode1d * nnode1d - 1;
466  break;
467 
468  default:
469 
470  std::ostringstream error_stream;
471  error_stream << "Never get here. vertex: " << Direct_string[vertex]
472  << std::endl;
473  throw OomphLibError(
474  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
475  break;
476  }
477 
478 
479  std::ostringstream error_stream;
480  error_stream << "Never get here. vertex: " << Direct_string[vertex]
481  << std::endl;
482  throw OomphLibError(
483  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
484  }

References Direct_string, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RUB, and oomph::OcTreeNames::RUF.

Referenced by oomph::OcTreeForest::construct_up_right_equivalents(), and oomph::HelmholtzMGPreconditioner< DIM >::maximum_edge_width().

Member Data Documentation

◆ Colour

Vector< std::string > oomph::OcTree::Colour
staticprivate

Colours for neighbours in various directions.

Referenced by doc_face_neighbours(), doc_true_edge_neighbours(), and setup_static_data().

◆ Common_face

DenseMatrix< int > oomph::OcTree::Common_face
staticprivate

Determine common face of edges or octants. Slightly bizarre lookup scheme from Samet's book.

Referenced by gteq_edge_neighbour(), and setup_static_data().

◆ Cosi

Vector< int > oomph::OcTree::Cosi
staticprivate

Entry in rotation matrix: cos(i*90)

Referenced by construct_rotation_matrix(), and setup_static_data().

◆ Direct_string

Vector< std::string > oomph::OcTree::Direct_string
static

Translate (enumerated) directions into strings.

Referenced by oomph::OcTreeRoot::add_edge_neighbour_pt(), doc_face_neighbours(), doc_true_edge_neighbours(), edge_neighbour_is_face_neighbour(), oomph::OcTreeRoot::edge_neighbour_pt(), faces_of_common_edge(), gteq_edge_neighbour(), gteq_face_neighbour(), gteq_true_edge_neighbour(), oomph::RefineableQElement< 3 >::interpolated_zeta_on_face(), oomph::EighthSphereDomain::macro_element_boundary(), oomph::ExtrudedDomain::macro_element_boundary(), oomph::OcTreeRoot::nedge_neighbour(), oomph::OcTreeForest::oc_face_neigh_pt(), rotate(), setup_static_data(), and vertex_to_node_number().

◆ Direction_to_vector

Vector< Vector< int > > oomph::OcTree::Direction_to_vector
static

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)

Vector storing the information to translate a direction into a vector of directions (in the three coordinate directions)

Referenced by oomph::PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >::further_build(), oomph::RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >::further_build(), oomph::PRefineableQCrouzeixRaviartElement< DIM >::further_build(), oomph::RefineableQPVDElementWithPressure< DIM >::further_build(), oomph::RefineableQElement< 3 >::get_bcs(), oomph::RefineableQElement< 3 >::get_boundaries(), oomph::RefineableSolidQElement< 3 >::get_solid_bcs(), get_the_other_face(), gteq_face_neighbour(), gteq_true_edge_neighbour(), rotate(), oomph::RefineableQElement< 3 >::setup_father_bounds(), and setup_static_data().

◆ Is_adjacent

DenseMatrix< bool > oomph::OcTree::Is_adjacent
staticprivate

Array of direction/octant adjacency scheme: Is_adjacent(direction,octant): Is face/edge direction adjacent to octant octant ? (Table in Samet's book)

Array of direction/octant adjacency scheme: Is_adjacent(direction,octant): Is face/edge direction adjacent to octant octant ? Table in Samet's book.

Referenced by gteq_edge_neighbour(), gteq_face_neighbour(), and setup_static_data().

◆ Reflect

DenseMatrix< int > oomph::OcTree::Reflect
staticprivate

Reflection scheme: Reflect(direction,octant): Get mirror of octant/edge in specified direction. E.g. Reflect(LDF,L)=RDF

Referenced by gteq_edge_neighbour(), gteq_face_neighbour(), and setup_static_data().

◆ Reflect_edge

Vector< int > oomph::OcTree::Reflect_edge
static

Get opposite edge, e.g. Reflect_edge[DB]=UF.

Referenced by gteq_true_edge_neighbour(), and setup_static_data().

◆ Reflect_face

Vector< int > oomph::OcTree::Reflect_face
static

Get opposite face, e.g. Reflect_face[L]=R.

Referenced by oomph::OcTreeForest::construct_up_right_equivalents(), gteq_face_neighbour(), and setup_static_data().

◆ Reflect_vertex

Vector< int > oomph::OcTree::Reflect_vertex
static

Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF.

Referenced by setup_static_data().

◆ S_base

DenseMatrix< double > oomph::OcTree::S_base
staticprivate

s_base(i,direction): Initial value for coordinate s[i] on the face indicated by direction (L/R/U/D/F/B)

Referenced by doc_face_neighbours(), gteq_face_neighbour(), and setup_static_data().

◆ S_base_edge

DenseMatrix< double > oomph::OcTree::S_base_edge
staticprivate

S_base_edge(i,edge): Initial value for coordinate s[i] on the specified edge (LF/RF/...).

Referenced by doc_true_edge_neighbours(), gteq_true_edge_neighbour(), and setup_static_data().

◆ S_direct_edge

DenseMatrix< double > oomph::OcTree::S_direct_edge
staticprivate

Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son specified by son_octant has an offset. If we project the son_octant's left/down/back vertex onto the father's edge edge, it is located at the in-face coordinate s_lo = h/2 S_direct_edge(edge,son_octant).

Referenced by gteq_edge_neighbour(), and setup_static_data().

◆ S_directhi

DenseMatrix< double > oomph::OcTree::S_directhi
staticprivate

Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son specified by son_octant has an offset. If we project the son_octant's left/down/back vertex onto the father's face face, it is located at the in-face coordinate s_hi = h/2 S_directlhi(face,son_octant). [See discussion of s_steplo for an explanation of the subscripts _hi and _lo.]

Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son specified by son_octant has an offset. If we project the son_octant's left/down/back vertex onto the father's face face, it is located at the in-face coordinate s_hi = h/2 S_directlhi(face,son_octant). [See discussion of S_steplo for an explanation of the subscripts _hi and _lo.]

Referenced by gteq_face_neighbour(), and setup_static_data().

◆ S_directlo

DenseMatrix< double > oomph::OcTree::S_directlo
staticprivate

Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son specified by son_octant has an offset. If we project the son_octant's left/down/back vertex onto the father's face face, it is located at the in-face coordinate s_lo = h/2 S_directlo(face,son_octant). [See discussion of s_steplo for an explanation of the subscripts _hi and _lo.]

Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son specified by son_octant has an offset. If we project the son_octant's left/down/back vertex onto the father's face face, it is located at the in-face coordinate s_lo = h/2 S_directlo(face,son_octant). [See discussion of S_steplo for an explanation of the subscripts _hi and _lo.]

Referenced by gteq_face_neighbour(), and setup_static_data().

◆ S_step_edge

DenseMatrix< double > oomph::OcTree::S_step_edge
staticprivate

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 edge [DB,UB,...], then an increase in s from -1 to +1 corresponds to a change of s_step_edge(i,edge) in the 3D coordinates s[i].

Referenced by doc_true_edge_neighbours(), gteq_true_edge_neighbour(), and setup_static_data().

◆ S_stephi

DenseMatrix< double > oomph::OcTree::S_stephi
staticprivate

If we're located on face face [L/R/F/B/U/D], then an increase in s_hi from -1 to +1 corresponds to a change of s_stephi(i,face) in the 3D coordinate \ s[i]. [Read the discussion of s_steplo for an explanation of the subscripts _hi and _lo.]

If we're located on face face [L/R/F/B/U/D], then an increase in s_hi from -1 to +1 corresponds to a change of S_stephi(i,face) in the 3D coordinate \ s[i]. [Read the discussion of S_steplo for an explanation of the subscripts _hi and _lo.]

Referenced by doc_face_neighbours(), gteq_face_neighbour(), and setup_static_data().

◆ S_steplo

DenseMatrix< double > oomph::OcTree::S_steplo
staticprivate

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 _lo and the one with the larger (3D) index with the subscript _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 face [L/R/F/B/U/D], then an increase in s_lo from -1 to +1 corresponds to a change of s_steplo(i,face) in the 3D coordinate s[i].

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 _lo and the one with the larger (3D) index with the subscript _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 face [L/R/F/B/U/D], then an increase in s_lo from -1 to +1 corresponds to a change of S_steplo(i,face) in the 3D coordinate s[i]. S_steplo(i,direction)

Referenced by doc_face_neighbours(), gteq_face_neighbour(), and setup_static_data().

◆ Sini

Vector< int > oomph::OcTree::Sini
staticprivate

Entry in rotation matrix sin(i*90)

Entry in rotation matrix: sin(i*90)

Referenced by construct_rotation_matrix(), and setup_static_data().

◆ Static_data_has_been_setup

bool oomph::OcTree::Static_data_has_been_setup = false
staticprotected

Bool indicating that static member data has been setup.

Referenced by oomph::OcTreeRoot::OcTreeRoot(), and setup_static_data().

◆ Up_and_right_equivalent_for_pairs_of_vertices

std::map< std::pair< std::pair< int, int >, std::pair< int, int > >, std::pair< int, int > > oomph::OcTree::Up_and_right_equivalent_for_pairs_of_vertices
static

Storage for the up/right-equivalents corresponding to two pairs of vertices along an element edge:

  • The first pair contains
    1. the vertex in the reference element
    2. 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
    1. the vertex at the other end of the edge in the reference element
    2. 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.

Referenced by oomph::OcTreeForest::construct_up_right_equivalents(), and setup_static_data().

◆ Vector_to_direction

std::map< Vector< int >, int > oomph::OcTree::Vector_to_direction
static

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

Map storing the information to translate a vector of directions (in the three coordinate directions) into a direction

Referenced by oomph::RefineableQElement< 3 >::get_bcs(), oomph::RefineableSolidQElement< 3 >::get_solid_bcs(), get_the_other_face(), rotate(), oomph::RefineableQElement< 3 >::setup_father_bounds(), and setup_static_data().

◆ Vertex_at_end_of_edge

Vector< Vector< int > > oomph::OcTree::Vertex_at_end_of_edge
static

Map of vectors containing the two vertices for each edge.

Vector of vectors containing the two vertices for each edge, e.g. Vertex_at_end_of_edge[LU][0]=LUB and Vertex_at_end_of_edge[LU][1]=LUF.

Referenced by oomph::OcTreeForest::construct_up_right_equivalents(), and setup_static_data().

The documentation for this class was generated from the following files: