octree.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #include <algorithm>
27 
28 #include "octree.h"
30 
31 namespace oomph
32 {
33  //========================================================================
34  /// Bool indicating that static member data has been setup
35  //========================================================================
37 
38 
39  // Public static member data:
40  //---------------------------
41 
42  //====================================================================
43  /// Translate (enumerated) directions into strings
44  //====================================================================
45  Vector<std::string> OcTree::Direct_string;
46 
47  //====================================================================
48  /// Get opposite face, e.g. Reflect_face[L]=R
49  //====================================================================
50  Vector<int> OcTree::Reflect_face;
51 
52  //====================================================================
53  /// Get opposite edge, e.g. Reflect_edge[DB]=UF
54  //====================================================================
55  Vector<int> OcTree::Reflect_edge;
56 
57  //====================================================================
58  /// Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF
59  //====================================================================
60  Vector<int> OcTree::Reflect_vertex;
61 
62  //====================================================================
63  /// Map of vectors containing the two vertices for each edge
64  //====================================================================
65  Vector<Vector<int>> OcTree::Vertex_at_end_of_edge;
66 
67  //====================================================================
68  /// Map storing the information to translate a vector of directions
69  /// (in the three coordinate directions) into a direction
70  //====================================================================
71  std::map<Vector<int>, int> OcTree::Vector_to_direction;
72 
73  //====================================================================
74  /// Vector storing the information to translate a direction into a
75  /// vector of directions (in the three coordinate directions)
76  //====================================================================
77  Vector<Vector<int>> OcTree::Direction_to_vector;
78 
79  //====================================================================
80  /// Storage for the up/right-equivalents corresponding to two
81  /// pairs of vertices along an element edge:
82  /// - The first pair contains
83  /// -# the vertex in the reference element
84  /// -# the corresponding vertex in the edge neighbour (i.e. the
85  /// vertex in the edge neighbour that is located at the same
86  /// position as that first vertex).
87  /// .
88  /// - The second pair contains
89  /// -# the vertex at the other end of the edge in the reference element
90  /// -# the corresponding vertex in the edge neighbour.
91  /// .
92  /// .
93  /// These two pairs completely define the relative rotation
94  /// between the reference element and its edge neighbour. The map
95  /// returns a pair which contains the up_equivalent and the
96  /// right_equivalent of the edge neighbour, i.e. it tells us
97  /// which direction in the edge neighbour coincides with the
98  /// up (or right) direction in the reference element.
99  //====================================================================
100  std::map<std::pair<std::pair<int, int>, std::pair<int, int>>,
101  std::pair<int, int>>
103 
104 
105  // Protected static member data:
106  //------------------------------
107 
108 
109  //====================================================================
110  /// Entry in rotation matrix: cos(i*90)
111  //====================================================================
112  Vector<int> OcTree::Cosi;
113 
114 
115  //====================================================================
116  /// Entry in rotation matrix: sin(i*90)
117  //====================================================================
118  Vector<int> OcTree::Sini;
119 
120  //====================================================================
121  /// Array of direction/octant adjacency scheme:
122  /// Is_adjacent(direction,octant): Is face/edge \c direction
123  /// adjacent to octant \c octant ? Table in Samet's book.
124  //====================================================================
125  DenseMatrix<bool> OcTree::Is_adjacent;
126 
127  //====================================================================
128  /// Reflection scheme: Reflect(direction,octant): Get mirror
129  /// of octant/edge in specified direction. E.g. Reflect(LDF,L)=RDF
130  //====================================================================
131  DenseMatrix<int> OcTree::Reflect;
132 
133  //====================================================================
134  /// Determine common face of edges or octants.
135  /// Slightly bizarre lookup scheme from Samet's book.
136  //====================================================================
137  DenseMatrix<int> OcTree::Common_face;
138 
139  //====================================================================
140  /// Colours for neighbours in various directions
141  //====================================================================
142  Vector<std::string> OcTree::Colour;
143 
144  //====================================================================
145  /// s_base(i,direction): Initial value for coordinate s[i] on
146  /// the face indicated by direction (L/R/U/D/F/B)
147  //====================================================================
148  DenseMatrix<double> OcTree::S_base;
149 
150  //====================================================================
151  /// Each face of the RefineableQElement<3> that is represented
152  /// by the octree is parametrised by two (of the three)
153  /// local coordinates that parametrise the entire 3D element. E.g.
154  /// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
155  /// is parametrised by (s[0],s[2]); etc. We always identify the
156  /// in-face coordinate with the lower (3D) index with the subscript
157  /// \c _lo and the one with the larger (3D) index with the subscript \c _hi.
158  /// Here we set up the translation scheme between the 2D in-face
159  /// coordinates (s_lo,s_hi) and the corresponding 3D coordinates:
160  /// If we're located on face \c face [L/R/F/B/U/D], then
161  /// an increase in s_lo from -1 to +1 corresponds to a change
162  /// of \c S_steplo(i,face) in the 3D coordinate \c s[i]. S_steplo(i,direction)
163  //====================================================================
164  DenseMatrix<double> OcTree::S_steplo;
165 
166 
167  //====================================================================
168  /// If we're located on face \c face [L/R/F/B/U/D], then
169  /// an increase in s_hi from -1 to +1 corresponds to a change
170  /// of \c S_stephi(i,face) in the 3D coordinate \ s[i].
171  /// [Read the discussion of \c S_steplo for an explanation of
172  /// the subscripts \c _hi and \c _lo.]
173  //====================================================================
174  DenseMatrix<double> OcTree::S_stephi;
175 
176  //====================================================================
177  /// Relative to the left/down/back vertex in any (father) octree, the
178  /// corresponding vertex in the son specified by \c son_octant has an offset.
179  /// If we project the son_octant's left/down/back vertex onto the
180  /// father's face \c face, it is located at the in-face coordinate
181  /// \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
182  /// \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
183  //====================================================================
184  DenseMatrix<double> OcTree::S_directlo;
185 
186  //====================================================================
187  /// Relative to the left/down/back vertex in any (father) octree, the
188  /// corresponding vertex in the son specified by \c son_octant has an offset.
189  /// If we project the son_octant's left/down/back vertex onto the
190  /// father's face \c face, it is located at the in-face coordinate
191  /// \c s_hi = h/2 \c S_directlhi(face,son_octant). [See discussion of
192  /// \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
193  //====================================================================
194  DenseMatrix<double> OcTree::S_directhi;
195 
196  //====================================================================
197  /// S_base_edge(i,edge): Initial value for coordinate s[i] on
198  /// the specified edge (LF/RF/...).
199  //====================================================================
200  DenseMatrix<double> OcTree::S_base_edge;
201 
202  //====================================================================
203  /// Each edge of the RefineableQElement<3> that is represented
204  /// by the octree is parametrised by one (of the three)
205  /// local coordinates that parametrise the entire 3D element.
206  /// If we're located on edge \c edge [DB,UB,...], then
207  /// an increase in s from -1 to +1 corresponds to a change
208  /// of \c s_step_edge(i,edge) in the 3D coordinates \c s[i].
209  //====================================================================
210  DenseMatrix<double> OcTree::S_step_edge;
211 
212  //====================================================================
213  /// Relative to the left/down/back vertex in any (father) octree, the
214  /// corresponding vertex in the son specified by \c son_octant has an offset.
215  /// If we project the son_octant's left/down/back vertex onto the
216  /// father's edge \c edge, it is located at the in-face coordinate
217  /// \c s_lo = h/2 \c S_direct_edge(edge,son_octant).
218  //====================================================================
219  DenseMatrix<double> OcTree::S_direct_edge;
220 
221 
222  // End static data
223  //----------------
224 
225 
226  //===================================================================
227  /// This function is used to translate the position of a vertex node
228  /// (given by his local number n into a vector giving the position of
229  /// this node in the local coordinates system.
230  /// It also needs the value of nnode1d to work.
231  //===================================================================
233  const unsigned& nnode1d)
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  }
264 
265  //==================================================================
266  /// Given an edge, this function returns the faces on which it lies.
267  //==================================================================
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
408 
409 
410  //==================================================================
411  /// Return the local node number of given vertex
412  /// in an element with nnode1d nodes in each coordinate direction.
413  //==================================================================
414  unsigned OcTree::vertex_to_node_number(const int& vertex,
415  const unsigned& nnode1d)
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  }
485 
486 
487  //==================================================================
488  /// Return the vertex of local (vertex) node n
489  /// in an element with nnode1d nodes in each coordinate direction.
490  //==================================================================
491  int OcTree::node_number_to_vertex(const unsigned& n, const unsigned& nnode1d)
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  }
555 
556 
557  //==================================================================
558  /// This function takes as argument two node numbers of two nodes
559  /// delimiting an edge, and one face of this edge and returns the
560  /// other face that is sharing this edge. The node numbers given to
561  /// this function MUST be vertices nodes to work.
562  /// it also need the value of nnode1d to work.
563  /// (\c face is a direction in the set U,D,F,B,L,R).
564  //==================================================================
565  int OcTree::get_the_other_face(const unsigned& n1,
566  const unsigned& n2,
567  const unsigned& nnode1d,
568  const int& face)
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  }
603 
604 
605  //====================================================================
606  /// Build the rotation matrix for a rotation around the axis \c axis of
607  /// an angle \c angle*90
608  //====================================================================
610  int& angle,
611  DenseMatrix<int>& mat)
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  }
660 
661  //===================================================================
662  /// Helper: Performs the operation Mat3=Mat1*Mat2
663  //===================================================================
665  const DenseMatrix<int>& mat2,
666  DenseMatrix<int>& mat3)
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  }
682 
683 
684  //===================================================================
685  /// Helper: Performs the operation Vect2=Mat*Vect1
686  //===================================================================
688  const Vector<int>& vect1,
689  Vector<int>& vect2)
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  }
702 
703 
704  //===================================================================
705  /// A rotation is defined by the newUp and newRight
706  /// directions; so if Up becomes newUp and Right becomes newRight
707  /// then dir becomes rotate(newUp,newRight,dir);
708  //===================================================================
709  int OcTree::rotate(const int& new_up, const int& new_right, const int& 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  }
744 
745 
746  //==================================================================
747  /// This function rotates a vector according to a rotation of
748  /// the axes that changes up to new_up and right to new_right.
749  //==================================================================
750  Vector<int> OcTree::rotate(const int& new_up,
751  const int& new_right,
752  const Vector<int>& dir)
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  }
1035 
1036 
1037  //====================================================================
1038  /// Setup static data for OcTree
1039  //====================================================================
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
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
2139 
2140  Vertex_at_end_of_edge[UB].resize(2);
2141  Vertex_at_end_of_edge[UB][0] = LUB; // Pattern: both other indices
2143 
2144 
2145  Vertex_at_end_of_edge[LB].resize(2);
2146  Vertex_at_end_of_edge[LB][0] = LUB; // Pattern: both other indices
2148 
2149  Vertex_at_end_of_edge[RB].resize(2);
2150  Vertex_at_end_of_edge[RB][0] = RUB; // Pattern: both other indices
2152 
2153 
2154  Vertex_at_end_of_edge[LD].resize(2);
2155  Vertex_at_end_of_edge[LD][0] = LDF; // Pattern: both other indices
2157 
2158  Vertex_at_end_of_edge[RD].resize(2);
2159  Vertex_at_end_of_edge[RD][0] = RDF; // Pattern: both other indices
2161 
2162  Vertex_at_end_of_edge[LU].resize(2);
2163  Vertex_at_end_of_edge[LU][0] = LUF; // Pattern: both other indices
2165 
2166  Vertex_at_end_of_edge[RU].resize(2);
2167  Vertex_at_end_of_edge[RU][0] = RUF; // Pattern: both other indices
2169 
2170 
2171  Vertex_at_end_of_edge[DF].resize(2);
2172  Vertex_at_end_of_edge[DF][0] = LDF; // Pattern: both other indices
2174 
2175  Vertex_at_end_of_edge[UF].resize(2);
2176  Vertex_at_end_of_edge[UF][0] = LUF; // Pattern: both other indices
2178 
2179  Vertex_at_end_of_edge[LF].resize(2);
2180  Vertex_at_end_of_edge[LF][0] = LUF; // Pattern: both other indices
2182 
2183  Vertex_at_end_of_edge[RF].resize(2);
2184  Vertex_at_end_of_edge[RF][0] = RUF; // Pattern: both other indices
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:
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
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  }
2763 
2764 
2765  //================================================================
2766  /// Is the edge neighbour (for edge "edge") specified via the pointer
2767  /// also a face neighbour for one of the two adjacent faces?
2768  //================================================================
2770  OcTree* edge_neigh_pt) const
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  }
3328 
3329 
3330  //================================================================
3331  /// Find (pointer to) `greater-or-equal-sized face neighbour' in
3332  /// given direction (L/R/U/D/F/B).
3333  /// Another way of interpreting this is that we're looking for
3334  /// the neighbour across the present element's face 'direction'.
3335  /// The various arguments return additional information about the
3336  /// size and relative orientation of the neighbouring octree.
3337  /// To interpret these we use the following
3338  /// <B>General convention:</B>
3339  /// - Each face of the element that is represented by the octree
3340  /// is parametrised by two (of the three)
3341  /// local coordinates that parametrise the entire 3D element. E.g.
3342  /// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
3343  /// is parametrised by (s[0],s[2]); etc. We always identify the
3344  /// in-face coordinate with the lower (3D) index with the subscript
3345  /// _lo and the one with the larger (3D) index with the subscript _hi.
3346  /// .
3347  /// With this convention, the interpretation of the arguments is
3348  /// as follows:
3349  /// - The vector \c translate_s turns the index of the local coordinate
3350  /// in the present octree into that of the neighbour. If there are no
3351  /// rotations then \c translate_s[i] = i.
3352  /// - In the present octree, the "south west" vertex of the face
3353  /// between the present octree and its neighbour is located at
3354  /// S_lo=-1, S_hi=-1. This point is located at the (3D) local
3355  /// coordinates (\c s_sw[0], \c s_sw[1], \c s_sw[2])
3356  /// in the neighbouring octree.
3357  /// - ditto with s_ne: In the present octree, the "north east" vertex
3358  /// of the face between the present octree and its neighbour is located at
3359  /// S_lo=+1, S_hi=+1. This point is located
3360  /// at the (3D) local coordinates (\c s_ne[0], \c s_ne[1], \c s_ne[2])
3361  /// in the neighbouring octree.
3362  /// - We're looking for a neighbour in the specified \c direction. When
3363  /// viewed from the neighbouring octree, the face that separates
3364  /// the present octree from its neighbour is the neighbour's face
3365  /// \c face. If there's no rotation between the two octrees, this is a
3366  /// simple reflection: For instance, if we're looking
3367  /// for a neighhbour in the \c R [ight] \c direction, \c face will
3368  /// be \c L [eft]
3369  /// - \c diff_level <= 0 indicates the difference in refinement levels between
3370  /// the two neighbours. If \c diff_level==0, the neighbour has the
3371  /// same size as the current octree.
3372  //=================================================================
3373  OcTree* OcTree::gteq_face_neighbour(const int& direction,
3374  Vector<unsigned>& translate_s,
3375  Vector<double>& s_sw,
3376  Vector<double>& s_ne,
3377  int& face,
3378  int& diff_level,
3379  bool& in_neighbouring_tree) const
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  }
3562 
3563  //================================================================
3564  /// Find (pointer to) `greater-or-equal-sized true edge neighbour' in
3565  /// the given direction (LB,RB,DB,UB [the back edges],
3566  /// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]).
3567  /// Another way of interpreting this is that we're looking for
3568  /// the neighbour across the present element's edge 'direction'.
3569  /// The various arguments return additional information about the
3570  /// size and relative orientation of the neighbouring octree.
3571  /// Each edge of the element that is represented by the octree
3572  /// is parametrised by one (of the three) local coordinates that
3573  /// parametrise the entire 3D element. E.g. the L[eft]B[ack] edge
3574  /// is parametrised by s[1]; the "low" vertex of this edge
3575  /// (located at the low value of this coordinate, i.e. at s[1]=-1)
3576  /// is L[eft]D[own]B[ack]. The "high" vertex of this edge (located
3577  /// at the high value of this coordinate, i.e. at s[1]=1) is
3578  /// L[eft]U[p]B[ack]; etc
3579  /// The interpretation of the arguments is as follows:
3580  /// - In a forest, an OcTree can have multiple edge neighbours
3581  /// (across an edge where multiple trees meet). \c i_root_edge_neighbour
3582  /// specifies which of these is used. Use this as "reverse communication":
3583  /// First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
3584  /// initialised to anything you want (zero, ideally). On return from
3585  /// the fct, \c n_root_edge_neighour contains the total number of true
3586  /// edge neighbours, so additional calls to the fct with
3587  /// \c i_root_edge_neighbour>0 can be made until they've all been visited.
3588  /// - The vector \c translate_s turns the index of the local coordinate
3589  /// in the present octree into that of the neighbour. If there are no
3590  /// rotations then \c translate_s[i] = i.
3591  /// - The "low" vertex of the edge in the present octree
3592  /// coincides with a certain vertex in the edge neighbour.
3593  /// In terms of the neighbour's local coordinates, this point is
3594  /// located at the (3D) local coordinates (\c s_lo[0], \c s_lo[1],
3595  /// \c s_lo[2])
3596  /// - ditto with s_hi: The "high" vertex of the edge in the present octree
3597  /// coincides with a certain vertex in the edge neighbour.
3598  /// In terms of the neighbour's local coordinates, this point is
3599  /// located at the (3D) local coordinates (\c s_hi[0], \c s_hi[1],
3600  /// \c s_hi[2])
3601  /// - We're looking for a neighbour in the specified \c direction. When
3602  /// viewed from the neighbouring octree, the edge that separates
3603  /// the present octree from its neighbour is the neighbour's edge
3604  /// \c edge. If there's no rotation between the two octrees, this is a
3605  /// simple reflection: For instance, if we're looking
3606  /// for a neighhbour in the \c DB \c direction, \c edge will
3607  /// be \c UF.
3608  /// - \c diff_level <= 0 indicates the difference in refinement levels between
3609  /// the two neighbours. If \c diff_level==0, the neighbour has the
3610  /// same size as the current octree.
3611  /// .
3612  /// \b Important: We're only looking for \b true edge neighbours
3613  /// i.e. edge neigbours that are not also face neighbours. This is an
3614  /// important difference to Samet's terminology. If the neighbour
3615  /// in a certain direction is not a true edge neighbour, or if there
3616  /// is no neighbour, then this function returns NULL.
3617  //=================================================================
3619  const int& direction,
3620  const unsigned& i_root_edge_neighbour,
3621  unsigned& nroot_edge_neighbour,
3622  Vector<unsigned>& translate_s,
3623  Vector<double>& s_lo,
3624  Vector<double>& s_hi,
3625  int& edge,
3626  int& diff_level) const
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  }
3809 
3810 
3811  //==========================================================================
3812  /// Find `greater-or-equal-sized face neighbour' in given direction
3813  /// (L/R/U/D/B/F).
3814  ///
3815  /// This is an auxiliary routine which allows neighbour finding in adjacent
3816  /// octrees. Needs to keep track of previous son types and
3817  /// the maximum level to which search is performed.
3818  ///
3819  /// Parameters:
3820  ///
3821  /// - direction: (L/R/U/D/B/F) Direction in which neighbour has to be found.
3822  /// - s_difflo/s_diffhi: Offset of left/down/back vertex from
3823  /// corresponding vertex in neighbour. Note that this is input/output
3824  /// as it needs to be incremented/decremented during the recursive calls
3825  /// to this function.
3826  /// - face: We're looking for the neighbour across our face 'direction'
3827  /// (L/R/U/D/B/F). When viewed from the neighbour, this face is
3828  /// `face' (L/R/U/D/B/F). [If there's no relative rotation between
3829  /// neighbours then this is a mere reflection, e.g. direction=F --> face=B
3830  /// etc.]
3831  /// - diff_level <= 0 indicates the difference in octree levels
3832  /// between the current element and its neighbour.
3833  /// - max_level is the maximum level to which the neighbour search is
3834  /// allowed to proceed. This is necessary because in a forest,
3835  /// the neighbour search isn't based on pure recursion.
3836  /// - orig_root_pt identifies the root node of the element whose
3837  /// neighbour we're really trying to find by all these recursive calls.
3838  //===========================================================================
3839  OcTree* OcTree::gteq_face_neighbour(const int& direction,
3840  double& s_difflo,
3841  double& s_diffhi,
3842  int& diff_level,
3843  bool& in_neighbouring_tree,
3844  int max_level,
3845  OcTreeRoot* orig_root_pt) const
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
3981 
3982 
3983  //==========================================================================
3984  /// Find `greater-or-equal-sized edge neighbour' in given direction
3985  /// (LB,RB,DB,UB [the back edges],
3986  /// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]).
3987  ///
3988  /// This is an auxiliary routine which allows neighbour finding in adjacent
3989  /// octrees. Needs to keep track of previous son types and
3990  /// the maximum level to which search is performed.
3991  ///
3992  /// Parameters:
3993  ///
3994  /// - direction: (LB,RB/...) Direction in which neighbour has to be found.
3995  /// - In a forest, an OcTree can have multiple edge neighbours
3996  /// (across an edge where multiple trees meet). \c i_root_edge_neighbour
3997  /// specifies which of these is used. Use this as "reverse communication":
3998  /// First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
3999  /// initialised to anything you want (zero, ideally). On return from
4000  /// the fct, \c n_root_edge_neighour contains the total number of true
4001  /// edge neighbours, so additional calls to the fct with
4002  /// \c i_root_edge_neighbour>0 can be made until they've all been visited.
4003  /// - s_diff: Offset of left/down/back vertex from
4004  /// corresponding vertex in
4005  /// neighbour. Note that this is input/output as it needs to be incremented/
4006  /// decremented during the recursive calls to this function.
4007  /// - diff_level <= 0 indicates the difference in octree levels
4008  /// between the current element and its neighbour.
4009  /// - max_level is the maximum level to which the neighbour search is
4010  /// allowed to proceed. This is necessary because in a forest,
4011  /// the neighbour search isn't based on pure recursion.
4012  /// - orig_root_pt identifies the root node of the element whose
4013  /// neighbour we're really trying to find by all these recursive calls.
4014  //===========================================================================
4015  OcTree* OcTree::gteq_edge_neighbour(const int& direction,
4016  const unsigned& i_root_edge_neighbour,
4017  unsigned& nroot_edge_neighbour,
4018  double& s_diff,
4019  int& diff_level,
4020  int max_level,
4021  OcTreeRoot* orig_root_pt) const
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
4190 
4191 
4192  //================================================================
4193  /// Self-test: Check neighbour finding routine. For each element
4194  /// in the tree and for each vertex, determine the
4195  /// distance between the vertex and its position in the
4196  /// neigbour. . If the difference is less than
4197  /// Tree::Max_neighbour_finding_tolerance.
4198  /// return success (0), otherwise failure (1)
4199  //=================================================================
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;
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  }
4263 
4264 
4265  //=================================================================
4266  /// Doc/check all face neighbours of octree (nodes) contained in the
4267  /// Vector forest_node_pt. Output into neighbours_file which can
4268  /// be viewed from tecplot with OcTreeNeighbours.mcr
4269  /// Neighbour info and errors are displayed on
4270  /// neighbours_txt_file. Finally, compute the max. error between
4271  /// vertices when viewed from neighbouring element.
4272  /// If the two filestreams are closed, output is suppressed.
4273  /// (Static function.)
4274  //=================================================================
4276  std::ofstream& neighbours_file,
4277  std::ofstream& neighbours_txt_file,
4278  double& max_error)
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  }
4536 
4537 
4538  /// ///////////////////////////////////////////////////////////////////
4539  /// ///////////////////////////////////////////////////////////////////
4540 
4541 
4542  //=================================================================
4543  /// Doc/check all true edge neighbours of octree (nodes) contained
4544  /// in the Vector forest_node_pt. Output into neighbours_file which can
4545  /// be viewed from tecplot with OcTreeNeighbours.mcr
4546  /// Neighbour info and errors are displayed on
4547  /// neighbours_txt_file. Finally, compute the max. error between
4548  /// vertices when viewed from neighbouring element.
4549  /// If the two filestreams are closed, output is suppressed.
4550  /// (Static function).
4551  //=================================================================
4553  std::ofstream& neighbours_file,
4554  std::ofstream& no_true_edge_file,
4555  std::ofstream& neighbours_txt_file,
4556  double& max_error)
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: &qu