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// Version 1.0; svn revision $LastChangedRevision: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #include<algorithm>
31 
32 #include "octree.h"
34 
35 
36 namespace oomph
37 {
38 
39 //========================================================================
40 /// Bool indicating that static member data has been setup
41 //========================================================================
43 
44 
45 
46 // Public static member data:
47 //---------------------------
48 
49 //====================================================================
50 /// Translate (enumerated) directions into strings
51 //====================================================================
52 Vector<std::string> OcTree::Direct_string;
53 
54 //====================================================================
55 /// Get opposite face, e.g. Reflect_face[L]=R
56 //====================================================================
57 Vector<int> OcTree::Reflect_face;
58 
59 //====================================================================
60 /// Get opposite edge, e.g. Reflect_edge[DB]=UF
61 //====================================================================
62 Vector<int> OcTree::Reflect_edge;
63 
64 //====================================================================
65 /// Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF
66 //====================================================================
67 Vector<int> OcTree::Reflect_vertex;
68 
69 //====================================================================
70 /// Map of vectors containing the two vertices for each edge
71 //====================================================================
72 Vector<Vector<int> > OcTree::Vertex_at_end_of_edge;
73 
74 //====================================================================
75 /// Map storing the information to translate a vector of directions
76 /// (in the three coordinate directions) into a direction
77 //====================================================================
78 std::map<Vector<int>, int > OcTree::Vector_to_direction;
79 
80 //====================================================================
81 /// Vector storing the information to translate a direction into a
82 /// vector of directions (in the three coordinate directions)
83 //====================================================================
84 Vector<Vector<int> > OcTree::Direction_to_vector;
85 
86 //====================================================================
87 /// \short Storage for the up/right-equivalents corresponding to two
88 /// pairs of vertices along an element edge:
89 /// - The first pair contains
90 /// -# the vertex in the reference element
91 /// -# the corresponding vertex in the edge neighbour (i.e. the
92 /// vertex in the edge neighbour that is located at the same
93 /// position as that first vertex).
94 /// .
95 /// - The second pair contains
96 /// -# the vertex at the other end of the edge in the reference element
97 /// -# the corresponding vertex in the edge neighbour.
98 /// .
99 /// .
100 /// These two pairs completely define the relative rotation
101 /// between the reference element and its edge neighbour. The map
102 /// returns a pair which contains the up_equivalent and the
103 /// right_equivalent of the edge neighbour, i.e. it tells us
104 /// which direction in the edge neighbour coincides with the
105 /// up (or right) direction in the reference element.
106 //====================================================================
107 std::map<std::pair<std::pair<int,int>,std::pair<int,int> >,std::pair<int,int> >
109 
110 
111 
112 
113 
114 // Protected static member data:
115 //------------------------------
116 
117 
118 //====================================================================
119  /// Entry in rotation matrix: cos(i*90)
120 //====================================================================
121 Vector<int> OcTree::Cosi;
122 
123 
124 //====================================================================
125  /// Entry in rotation matrix: sin(i*90)
126 //====================================================================
127 Vector<int> OcTree::Sini;
128 
129 //====================================================================
130 /// Array of direction/octant adjacency scheme:
131 /// Is_adjacent(direction,octant): Is face/edge \c direction
132 /// adjacent to octant \c octant ? Table in Samet's book.
133 //====================================================================
134 DenseMatrix<bool> OcTree::Is_adjacent;
135 
136 //====================================================================
137 /// Reflection scheme: Reflect(direction,octant): Get mirror
138 /// of octant/edge in specified direction. E.g. Reflect(LDF,L)=RDF
139 //====================================================================
140 DenseMatrix<int> OcTree::Reflect;
141 
142 //====================================================================
143 /// Determine common face of edges or octants.
144 /// Slightly bizarre lookup scheme from Samet's book.
145 //====================================================================
146 DenseMatrix<int> OcTree::Common_face;
147 
148 //====================================================================
149 /// Colours for neighbours in various directions
150 //====================================================================
151 Vector<std::string> OcTree::Colour;
152 
153 //====================================================================
154 /// s_base(i,direction): Initial value for coordinate s[i] on
155 /// the face indicated by direction (L/R/U/D/F/B)
156 //====================================================================
157 DenseMatrix<double> OcTree::S_base;
158 
159 //====================================================================
160 /// Each face of the RefineableQElement<3> that is represented
161 /// by the octree is parametrised by two (of the three)
162 /// local coordinates that parametrise the entire 3D element. E.g.
163 /// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
164 /// is parametrised by (s[0],s[2]); etc. We always identify the
165 /// in-face coordinate with the lower (3D) index with the subscript
166 /// \c _lo and the one with the larger (3D) index with the subscript \c _hi.
167 /// Here we set up the translation scheme between the 2D in-face
168 /// coordinates (s_lo,s_hi) and the corresponding 3D coordinates:
169 /// If we're located on face \c face [L/R/F/B/U/D], then
170 /// an increase in s_lo from -1 to +1 corresponds to a change
171 /// of \c S_steplo(i,face) in the 3D coordinate \c s[i]. S_steplo(i,direction)
172 //====================================================================
173 DenseMatrix<double> OcTree::S_steplo;
174 
175 
176 //====================================================================
177 /// If we're located on face \c face [L/R/F/B/U/D], then
178 /// an increase in s_hi from -1 to +1 corresponds to a change
179 /// of \c S_stephi(i,face) in the 3D coordinate \ s[i].
180 /// [Read the discussion of \c S_steplo for an explanation of
181 /// the subscripts \c _hi and \c _lo.]
182 //====================================================================
183 DenseMatrix<double> OcTree::S_stephi;
184 
185 //====================================================================
186 /// Relative to the left/down/back vertex in any (father) octree, the
187 /// corresponding vertex in the son specified by \c son_octant has an offset.
188 /// If we project the son_octant's left/down/back vertex onto the
189 /// father's face \c face, it is located at the in-face coordinate
190 /// \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
191 /// \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
192 //====================================================================
193 DenseMatrix<double> OcTree::S_directlo;
194 
195 //====================================================================
196 /// Relative to the left/down/back vertex in any (father) octree, the
197 /// corresponding vertex in the son specified by \c son_octant has an offset.
198 /// If we project the son_octant's left/down/back vertex onto the
199 /// father's face \c face, it is located at the in-face coordinate
200 /// \c s_hi = h/2 \c S_directlhi(face,son_octant). [See discussion of
201 /// \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
202 //====================================================================
203 DenseMatrix<double> OcTree::S_directhi;
204 
205 //====================================================================
206 /// S_base_edge(i,edge): Initial value for coordinate s[i] on
207 /// the specified edge (LF/RF/...).
208 //====================================================================
209 DenseMatrix<double> OcTree::S_base_edge;
210 
211 //====================================================================
212 /// Each edge of the RefineableQElement<3> that is represented
213  /// by the octree is parametrised by one (of the three)
214  /// local coordinates that parametrise the entire 3D element.
215  /// If we're located on edge \c edge [DB,UB,...], then
216  /// an increase in s from -1 to +1 corresponds to a change
217  /// of \c s_step_edge(i,edge) in the 3D coordinates \c s[i].
218 //====================================================================
219 DenseMatrix<double> OcTree::S_step_edge;
220 
221 //====================================================================
222 /// Relative to the left/down/back vertex in any (father) octree, the
223 /// corresponding vertex in the son specified by \c son_octant has an offset.
224 /// If we project the son_octant's left/down/back vertex onto the
225 /// father's edge \c edge, it is located at the in-face coordinate
226 /// \c s_lo = h/2 \c S_direct_edge(edge,son_octant).
227 //====================================================================
228 DenseMatrix<double> OcTree::S_direct_edge;
229 
230 
231 
232 
233 
234 // End static data
235 //----------------
236 
237 
238 
239 //===================================================================
240 /// This function is used to translate the position of a vertex node
241 /// (given by his local number n into a vector giving the position of
242 /// this node in the local coordinates system.
243 /// It also needs the value of nnode1d to work.
244 //===================================================================
246  const unsigned& nnode1d)
247 {
248 #ifdef PARANOID
249  if ((n!=0)&&
250  (n!=nnode1d-1)&&
251  (n!=(nnode1d-1)*nnode1d)&&
252  (n!=nnode1d*nnode1d-1)&&
253  (n!=(nnode1d*nnode1d)*(nnode1d-1)+0)&&
254  (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d-1)&&
255  (n!=(nnode1d*nnode1d)*(nnode1d-1)+(nnode1d-1)*nnode1d)&&
256  (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d*nnode1d-1))
257  {
258  std::ostringstream error_stream;
259  error_stream << "Node " << n
260  << " is not a vertex node in a brick element with "
261  << nnode1d << " nodes along each edge!";
262 
263  throw OomphLibError(error_stream.str(),
264  OOMPH_CURRENT_FUNCTION,
265  OOMPH_EXCEPTION_LOCATION);
266  }
267 #endif
268  int a,b,c;
269  Vector<int> result_vect(3);
270  a=n / (nnode1d*nnode1d);
271  b=(n-a*nnode1d*nnode1d)/nnode1d;
272  c=n-a*nnode1d*nnode1d-b*nnode1d;
273 
274  result_vect[0]=2*c/(nnode1d-1)-1;
275  result_vect[1]=2*b/(nnode1d-1)-1;
276  result_vect[2]=2*a/(nnode1d-1)-1;
277 
278  return result_vect;
279 }
280 
281 
282 
283 
284 
285 //==================================================================
286 /// Return the local node number of given vertex
287 /// in an element with nnode1d nodes in each coordinate direction.
288 //==================================================================
289 unsigned OcTree::vertex_to_node_number(const int& vertex,
290  const unsigned& nnode1d)
291 {
292  using namespace OcTreeNames;
293 
294 #ifdef PARANOID
295  if ((vertex!=LDB)&&(vertex!=RDB)&&
296  (vertex!=LUB)&&(vertex!=RUB)&&
297  (vertex!=LDF)&&(vertex!=RDF)&&
298  (vertex!=LUF)&&(vertex!=RUF))
299  {
300  std::ostringstream error_stream;
301  error_stream << "Wrong vertex: " << Direct_string[vertex]
302  << std::endl;
303  throw OomphLibError(error_stream.str(),
304  OOMPH_CURRENT_FUNCTION,
305  OOMPH_EXCEPTION_LOCATION);
306  }
307 #endif
308 
309  switch (vertex)
310  {
311 
312  case LDB:
313  return 0;
314  break;
315 
316  case RDB:
317  return nnode1d-1;
318  break;
319 
320 
321  case LUB:
322  return nnode1d*(nnode1d-1);
323  break;
324 
325  case RUB:
326  return nnode1d*nnode1d-1;
327  break;
328 
329 
330  case LDF:
331  return nnode1d*nnode1d*(nnode1d-1);
332  break;
333 
334 
335  case RDF:
336  return (nnode1d*nnode1d+1)*(nnode1d-1);
337  break;
338 
339  case LUF:
340  return nnode1d*nnode1d*nnode1d-nnode1d;
341  break;
342 
343  case RUF:
344  return nnode1d*nnode1d*nnode1d-1;
345  break;
346 
347  default:
348 
349  std::ostringstream error_stream;
350  error_stream << "Never get here. vertex: " << Direct_string[vertex]
351  << std::endl;
352  throw OomphLibError(error_stream.str(),
353  OOMPH_CURRENT_FUNCTION,
354  OOMPH_EXCEPTION_LOCATION);
355  break;
356  }
357 
358 
359 
360  std::ostringstream error_stream;
361  error_stream << "Never get here. vertex: " << Direct_string[vertex]
362  << std::endl;
363  throw OomphLibError(error_stream.str(),
364  OOMPH_CURRENT_FUNCTION,
365  OOMPH_EXCEPTION_LOCATION);
366 
367 
368 }
369 
370 
371 //==================================================================
372 /// \short Return the vertex of local (vertex) node n
373 /// in an element with nnode1d nodes in each coordinate direction.
374 //==================================================================
375  int OcTree::node_number_to_vertex(const unsigned& n,
376  const unsigned& nnode1d)
377 {
378  using namespace OcTreeNames;
379 
380 #ifdef PARANOID
381  if ((n!=0)&&
382  (n!=nnode1d-1)&&
383  (n!=(nnode1d-1)*nnode1d)&&
384  (n!=nnode1d*nnode1d-1)&&
385  (n!=(nnode1d*nnode1d)*(nnode1d-1)+0)&&
386  (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d-1)&&
387  (n!=(nnode1d*nnode1d)*(nnode1d-1)+(nnode1d-1)*nnode1d)&&
388  (n!=(nnode1d*nnode1d)*(nnode1d-1)+nnode1d*nnode1d-1))
389  {
390  std::ostringstream error_stream;
391  error_stream << "Node " << n
392  << " is not a vertex node in a brick element with "
393  << nnode1d << " nodes along each edge!";
394 
395  throw OomphLibError(error_stream.str(),
396  OOMPH_CURRENT_FUNCTION,
397  OOMPH_EXCEPTION_LOCATION);
398  }
399 #endif
400 
401  if (n==0) {return LDB;}
402  else if (n==nnode1d-1) {return RDB;}
403  else if (n==nnode1d*(nnode1d-1)) {return LUB;}
404  else if (n==nnode1d*nnode1d-1) {return RUB;}
405  else if (n==nnode1d*nnode1d*(nnode1d-1)) {return LDF;}
406  else if (n==(nnode1d*nnode1d+1)*(nnode1d-1)) {return RDF;}
407  else if (n==nnode1d*nnode1d*nnode1d-nnode1d) {return LUF;}
408  else if (n==nnode1d*nnode1d*nnode1d-1) {return RUF;}
409  else
410  {
411  std::ostringstream error_stream;
412  error_stream << "Never get here. local node number: " << n
413  << " is not a vertex node in a brick element with "
414  << nnode1d << " nodes along each edge!"
415  << std::endl;
416  throw OomphLibError(error_stream.str(),
417  OOMPH_CURRENT_FUNCTION,
418  OOMPH_EXCEPTION_LOCATION);
419  }
420 }
421 
422 
423 
424 
425 //==================================================================
426 /// This function takes as argument two node numbers of two nodes
427 /// delimiting an edge, and one face of this edge and returns the
428 /// other face that is sharing this edge. The node numbers given to
429 /// this function MUST be vertices nodes to work.
430 /// it also need the value of nnode1d to work.
431 /// (\c face is a direction in the set U,D,F,B,L,R).
432 //==================================================================
433 int OcTree::get_the_other_face(const unsigned& n1,
434  const unsigned& n2,
435  const unsigned& nnode1d,
436  const int& face)
437 {
438  Vector<int> vect_node1(3);
439  Vector<int> vect_node2(3);
440  Vector<int> vect_face(3);
441  Vector<int> vect_other_face(3);
442 
443  // translate the nodes to vectors
444  vect_node1=vertex_node_to_vector(n1,nnode1d);
445  vect_node2=vertex_node_to_vector(n2,nnode1d);
446  //translate the face to a vector
447  vect_face=Direction_to_vector[face];
448 
449  // compute the vector to the other face -- magic, courtesy of Renaud
450  // Schleck!
451  for(unsigned i=0;i<3;i++)
452  {
453  vect_other_face[i]=(vect_node1[i]+vect_node2[i])/2-vect_face[i];
454 
455 #ifdef PARANOID
456  if((vect_other_face[i]!=1)&&(vect_other_face[i]!=-1)&&
457  (vect_other_face[i]!=0))
458  {
459  throw OomphLibError("The nodes given are not vertices",
460  OOMPH_CURRENT_FUNCTION,
461  OOMPH_EXCEPTION_LOCATION);
462  }
463 #endif
464  }
465 
466  //return the corresponding face
467  return Vector_to_direction[vect_other_face];
468 }
469 
470 
471 
472 //====================================================================
473 /// Build the rotation matrix for a rotation around the axis \c axis of
474 /// an angle \c angle*90
475 //====================================================================
477  DenseMatrix<int>& mat)
478 {
479  using namespace OcTreeNames;
480  int a=0,b=0,c=0,i,j;
481 
482  switch(axis)
483  {
484  case R : a=1; b=2; c=0; break;
485  case U : a=2; b=0; c=1; break;
486  case F : a=0; b=1; c=2; break;
487  default :
488  throw OomphLibError("Bad Axis",
489  OOMPH_CURRENT_FUNCTION,
490  OOMPH_EXCEPTION_LOCATION);
491  }
492  for(i=0;i<3;i++)
493  {
494  for(j=0;j<3;j++)
495  {
496  mat(i,j)=0;
497  }
498  }
499  mat(a,a) = Cosi[angle];
500  mat(b,b) = Cosi[angle];
501  mat(a,b) = -Sini[angle];
502  mat(b,a) = Sini[angle];
503  mat(c,c) = 1;
504 }
505 
506 //===================================================================
507 /// Helper: Performs the operation Mat3=Mat1*Mat2
508 //===================================================================
510  const DenseMatrix<int>& mat2,
511  DenseMatrix<int>& mat3)
512 {
513  int Sum,i,j,k;
514  for(i=0;i<3;i++)
515  {
516  for(j=0;j<3;j++)
517  {
518  Sum=0;
519  for(k=0;k<3;k++)
520  {
521  Sum+=mat1(i,k)*mat2(k,j);
522  }
523  mat3(i,j)=Sum;
524  }
525  }
526 }
527 
528 
529 
530 //===================================================================
531 /// Helper: Performs the operation Vect2=Mat*Vect1
532 //===================================================================
534  const Vector<int>& vect1,
535  Vector<int>& vect2)
536 {
537  int i,k,sum;
538  for(i=0;i<3;i++)
539  {
540  sum=0;
541  for(k=0;k<3;k++)
542  {
543  sum+=mat(i,k)*vect1[k];
544  }
545  vect2[i]=sum;
546  }
547 }
548 
549 
550 
551 
552 //===================================================================
553 /// A rotation is defined by the newUp and newRight
554 /// directions; so if Up becomes newUp and Right becomes newRight
555 /// then dir becomes rotate(newUp,newRight,dir);
556 //===================================================================
557 int OcTree::rotate(const int& new_up, const int& new_right,
558  const int& dir)
559 {
560  using namespace OcTreeNames;
561 
562 #ifdef PARANOID
563  if ((new_up!=L)&&(new_up!=R)&&
564  (new_up!=F)&&(new_up!=B)&&
565  (new_up!=U)&&(new_up!=D))
566  {
567  std::ostringstream error_stream;
568  error_stream << "Wrong new_up: " << Direct_string[new_up]
569  << std::endl;
570  throw OomphLibError(error_stream.str(),
571  OOMPH_CURRENT_FUNCTION,
572  OOMPH_EXCEPTION_LOCATION);
573  }
574  if ((new_right!=L)&&(new_right!=R)&&
575  (new_right!=F)&&(new_right!=B)&&
576  (new_right!=U)&&(new_right!=D))
577  {
578  std::ostringstream error_stream;
579  error_stream << "Wrong new_right: " << Direct_string[new_right]
580  << std::endl;
581  throw OomphLibError(error_stream.str(),
582  OOMPH_CURRENT_FUNCTION,
583  OOMPH_EXCEPTION_LOCATION);
584  }
585 #endif
586 
587  Vector<int> vect_dir(3);
588  Vector<int> vect_new_dir(3);
589  // translate the direction to a vector
590  vect_dir=Direction_to_vector[dir];
591 
592  //rotate it
593  vect_new_dir=rotate(new_up,new_right,vect_dir);
594 
595  // translate the vector back to a direction
596  return Vector_to_direction[vect_new_dir];
597 }
598 
599 
600 //==================================================================
601 /// This function rotates a vector according to a rotation of
602 /// the axes that changes up to new_up and right to new_right.
603 //==================================================================
604 Vector<int> OcTree::rotate(const int& new_up,const int& new_right,
605  const Vector<int>& dir)
606 {
607  using namespace OcTreeNames;
608 
609 #ifdef PARANOID
610  if ((new_up!=L)&&(new_up!=R)&&
611  (new_up!=F)&&(new_up!=B)&&
612  (new_up!=U)&&(new_up!=D))
613  {
614  std::ostringstream error_stream;
615  error_stream << "Wrong new_up: " << Direct_string[new_up]
616  << std::endl;
617  throw OomphLibError(error_stream.str(),
618  OOMPH_CURRENT_FUNCTION,
619  OOMPH_EXCEPTION_LOCATION);
620  }
621  if ((new_right!=L)&&(new_right!=R)&&
622  (new_right!=F)&&(new_right!=B)&&
623  (new_right!=U)&&(new_right!=D))
624  {
625  std::ostringstream error_stream;
626  error_stream << "Wrong new_right: " << Direct_string[new_right]
627  << std::endl;
628  throw OomphLibError(error_stream.str(),
629  OOMPH_CURRENT_FUNCTION,
630  OOMPH_EXCEPTION_LOCATION);
631  }
632 #endif
633 
634 
635  // Every possible rotation of an element can be produced by the
636  // compostition of two (or one) rotations about one of the main
637  // axes of the element (R, U, F). So for each (newRight,newUp)
638  // we define the (angle1,axis1) of the first rotation and the
639  // (angle2, axis2) of the second one (if needed)
640 
641  int axis1,axis2,angle1,angle2,nrot;
642  nrot=2;
643  if (new_up == U)
644  {
645  nrot=1;
646  axis1=U;
647  switch(new_right)
648  {
649  case R : angle1=0; break;
650  case B : angle1=1; break;
651  case L : angle1=2; break;
652  case F : angle1=3; break;
653  default :
654  std::ostringstream error_stream;
655  error_stream << "new_right is "
656  << new_right << " (" << Direct_string[new_right]
657  << ") "
658  << "it should be R, B, L, or F" << std::endl;
659  throw OomphLibError(error_stream.str(),
660  OOMPH_CURRENT_FUNCTION,
661  OOMPH_EXCEPTION_LOCATION);
662  }
663  }
664 
665  if (new_up == D)
666  {
667  switch(new_right)
668  {
669  case R :
670  nrot=1;
671  axis1=R; angle1=2;
672  break;
673  case B :
674  axis1=R; angle1=2;
675  axis2=U; angle2=1;
676  break;
677  case L:
678  nrot=1;
679  axis1=F; angle1=2;
680  break;
681  case F :
682  axis1=R; angle1=2;
683  axis2=U; angle2=3;
684  break;
685  default :
686  std::ostringstream error_stream;
687  error_stream << "new_right is "
688  << new_right << " (" << Direct_string[new_right]
689  << ") "
690  << "it should be R, B, L, or F" << std::endl;
691  throw OomphLibError(error_stream.str(),
692  OOMPH_CURRENT_FUNCTION,
693  OOMPH_EXCEPTION_LOCATION);
694 
695  }
696  }
697 
698  if (new_up == R)
699  {
700  switch(new_right)
701  {
702  case D :
703  nrot=1;
704  axis1=F; angle1=3;
705  break;
706  case B :
707  axis1=F; angle1=3;
708  axis2=R; angle2=1;
709  break;
710  case U:
711  axis1=F; angle1=1;
712  axis2=U; angle2=2;
713  break;
714  case F:
715  axis1=F; angle1=3;
716  axis2=R; angle2=3;
717  break;
718  default :
719  std::ostringstream error_stream;
720  error_stream << "new_right is "
721  << new_right << " (" << Direct_string[new_right]
722  << ") "
723  << "it should be D, B, U, or F" << std::endl;
724  throw OomphLibError(error_stream.str(),
725  OOMPH_CURRENT_FUNCTION,
726  OOMPH_EXCEPTION_LOCATION);
727  }
728  }
729 
730  if (new_up == L)
731  {
732  switch(new_right)
733  {
734  case D :
735  axis1=F; angle1=1;
736  axis2=R; angle2=2;
737  break;
738  case B :
739  axis1=F; angle1=1;
740  axis2=R; angle2=3;
741  break;
742  case U :
743  nrot=1;
744  axis1=F; angle1=1;
745  break;
746  case F :
747  axis1=F; angle1=1;
748  axis2=R; angle2=1;
749  break;
750  default :
751  std::ostringstream error_stream;
752  error_stream << "new_right is "
753  << new_right << " (" << Direct_string[new_right]
754  << ") "
755  << "it should be D, B, U, or F" << std::endl;
756  throw OomphLibError(error_stream.str(),
757  OOMPH_CURRENT_FUNCTION,
758  OOMPH_EXCEPTION_LOCATION);
759  }
760  }
761 
762  if (new_up == F)
763  {
764  switch(new_right)
765  {
766  case R :
767  nrot=1;
768  axis1=R; angle1=1;
769  break;
770  case L :
771  axis1=R; angle1=1;
772  axis2=F; angle2=2;
773  break;
774  case U:
775  axis1=R; angle1=1;
776  axis2=F; angle2=1;
777  break;
778  case D:
779  axis1=R; angle1=1;
780  axis2=F; angle2=3;
781  break;
782  default :
783  std::ostringstream error_stream;
784  error_stream << "new_right is "
785  << new_right << " (" << Direct_string[new_right]
786  << ") "
787  << "it should be R, L, U, or D" << std::endl;
788  throw OomphLibError(error_stream.str(),
789  OOMPH_CURRENT_FUNCTION,
790  OOMPH_EXCEPTION_LOCATION);
791  }
792  }
793  if (new_up == B)
794  {
795  switch(new_right)
796  {
797  case R :
798  nrot=1;
799  axis1=R; angle1=3;
800  break;
801  case L :
802  axis1=R; angle1=3;
803  axis2=F; angle2=2;
804  break;
805  case U:
806  axis1=R; angle1=3;
807  axis2=F; angle2=1;
808  break;
809  case D:
810  axis1=R; angle1=3;
811  axis2=F; angle2=3;
812  break;
813  default :
814  std::ostringstream error_stream;
815  error_stream << "new_right is "
816  << new_right << " (" << Direct_string[new_right]
817  << ") "
818  << "it should be R, L, U, or D" << std::endl;
819  throw OomphLibError(error_stream.str(),
820  OOMPH_CURRENT_FUNCTION,
821  OOMPH_EXCEPTION_LOCATION);
822  }
823  }
824 
825  Vector<int> vect_new_dir(3);
826  DenseMatrix<int> mat1(3);
827  DenseMatrix<int> mat2(3);
828  DenseMatrix<int> mat3(3);
829 
830  // Then we build the first rotation matrix
831  construct_rotation_matrix(axis1,angle1,mat1);
832 
833  // If needed we build the second one
834  if (nrot==2)
835  {
836  // And we make the composition of the two rotations
837  // which is stored in Mat3
838  construct_rotation_matrix(axis2,angle2,mat2);
839  mult_mat_mat(mat2,mat1,mat3);
840  }
841  else
842  {
843  // Else we just copy Mat1 into Mat3
844  for(int i=0;i<3;i++)
845  {
846  for(int j=0;j<3;j++)
847  {
848  mat3(i,j)=mat1(i,j);
849  }
850  }
851  }
852 
853  // Rotate
854  mult_mat_vect(mat3,dir,vect_new_dir);
855 
856  // Return the corresponding vector
857  return vect_new_dir;
858 }
859 
860 
861 
862 
863 //====================================================================
864 /// Setup static data for OcTree
865 //====================================================================
867 {
868  using namespace OcTreeNames;
869 
870 
871 #ifdef PARANOID
873  {
874  std::ostringstream error_stream;
875  error_stream
876  << "Inconsistent enumeration! \n Tree::OMEGA=" << Tree::OMEGA
877  << "\nOcTree::OMEGA=" << OcTree::OMEGA
878  << std::endl;
879  throw OomphLibError(error_stream.str(),
880  OOMPH_CURRENT_FUNCTION,
881  OOMPH_EXCEPTION_LOCATION);
882  }
883 #endif
884 
885 // Private static data
886 //--------------------
887 
888  // Set flag to indicate that static data has been setup
890 
891  // Initialisation of cosi and sini used in rotation matrix
892  Cosi.resize(4);
893  Sini.resize(4);
894 
895  Cosi[0]=1;
896  Cosi[1]=0;
897  Cosi[2]=-1;
898  Cosi[3]=0;
899 
900  Sini[0]=0;
901  Sini[1]=1;
902  Sini[2]=0;
903  Sini[3]=-1;
904 
905 
906  // Build direction/octant adjacency scheme
907  // Is_adjacent(i_face,j_octant):
908  // Is face adjacent to octant?
909  // See the table in Samet book
910  Is_adjacent.resize(27,27);
911  Is_adjacent(L,LDB)=true;
912  Is_adjacent(R,LDB)=false;
913  Is_adjacent(D,LDB)=true;
914  Is_adjacent(U,LDB)=false;
915  Is_adjacent(B,LDB)=true;
916  Is_adjacent(F,LDB)=false;
917 
918  Is_adjacent(L,LDF)=true;
919  Is_adjacent(R,LDF)=false;
920  Is_adjacent(D,LDF)=true;
921  Is_adjacent(U,LDF)=false;
922  Is_adjacent(B,LDF)=false;
923  Is_adjacent(F,LDF)=true;
924 
925  Is_adjacent(L,LUB)=true;
926  Is_adjacent(R,LUB)=false;
927  Is_adjacent(D,LUB)=false;
928  Is_adjacent(U,LUB)=true;
929  Is_adjacent(B,LUB)=true;
930  Is_adjacent(F,LUB)=false;
931 
932  Is_adjacent(L,LUF)=true;
933  Is_adjacent(R,LUF)=false;
934  Is_adjacent(D,LUF)=false;
935  Is_adjacent(U,LUF)=true;
936  Is_adjacent(B,LUF)=false;
937  Is_adjacent(F,LUF)=true;
938 
939  Is_adjacent(L,RDB)=false;
940  Is_adjacent(R,RDB)=true;
941  Is_adjacent(D,RDB)=true;
942  Is_adjacent(U,RDB)=false;
943  Is_adjacent(B,RDB)=true;
944  Is_adjacent(F,RDB)=false;
945 
946  Is_adjacent(L,RDF)=false;
947  Is_adjacent(R,RDF)=true;
948  Is_adjacent(D,RDF)=true;
949  Is_adjacent(U,RDF)=false;
950  Is_adjacent(B,RDF)=false;
951  Is_adjacent(F,RDF)=true;
952 
953  Is_adjacent(L,RUB)=false;
954  Is_adjacent(R,RUB)=true;
955  Is_adjacent(D,RUB)=false;
956  Is_adjacent(U,RUB)=true;
957  Is_adjacent(B,RUB)=true;
958  Is_adjacent(F,RUB)=false;
959 
960  Is_adjacent(L,RUF)=false;
961  Is_adjacent(R,RUF)=true;
962  Is_adjacent(D,RUF)=false;
963  Is_adjacent(U,RUF)=true;
964  Is_adjacent(B,RUF)=false;
965  Is_adjacent(F,RUF)=true;
966 
967 
968  // Build direction/octant adjacency scheme
969  // Is_adjacent(i_edge,j_octant):
970  // Is edge adjacent to octant?
971  // See the table in Samet book
972  Is_adjacent(LD,LDB)=true;
973  Is_adjacent(LD,LDF)=true;
974  Is_adjacent(LD,LUB)=false;
975  Is_adjacent(LD,LUF)=false;
976  Is_adjacent(LD,RDB)=false;
977  Is_adjacent(LD,RDF)=false;
978  Is_adjacent(LD,RUB)=false;
979  Is_adjacent(LD,RUF)=false;
980 
981  Is_adjacent(LU,LDB)=false;
982  Is_adjacent(LU,LDF)=false;
983  Is_adjacent(LU,LUB)=true;
984  Is_adjacent(LU,LUF)=true;
985  Is_adjacent(LU,RDB)=false;
986  Is_adjacent(LU,RDF)=false;
987  Is_adjacent(LU,RUB)=false;
988  Is_adjacent(LU,RUF)=false;
989 
990 
991  Is_adjacent(LB,LDB)=true;
992  Is_adjacent(LB,LDF)=false;
993  Is_adjacent(LB,LUB)=true;
994  Is_adjacent(LB,LUF)=false;
995  Is_adjacent(LB,RDB)=false;
996  Is_adjacent(LB,RDF)=false;
997  Is_adjacent(LB,RUB)=false;
998  Is_adjacent(LB,RUF)=false;
999 
1000 
1001  Is_adjacent(LF,LDB)=false;
1002  Is_adjacent(LF,LDF)=true;
1003  Is_adjacent(LF,LUB)=false;
1004  Is_adjacent(LF,LUF)=true;
1005  Is_adjacent(LF,RDB)=false;
1006  Is_adjacent(LF,RDF)=false;
1007  Is_adjacent(LF,RUB)=false;
1008  Is_adjacent(LF,RUF)=false;
1009 
1010 
1011  Is_adjacent(RD,LDB)=false;
1012  Is_adjacent(RD,LDF)=false;
1013  Is_adjacent(RD,LUB)=false;
1014  Is_adjacent(RD,LUF)=false;
1015  Is_adjacent(RD,RDB)=true;
1016  Is_adjacent(RD,RDF)=true;
1017  Is_adjacent(RD,RUB)=false;
1018  Is_adjacent(RD,RUF)=false;
1019 
1020 
1021  Is_adjacent(RU,LDB)=false;
1022  Is_adjacent(RU,LDF)=false;
1023  Is_adjacent(RU,LUB)=false;
1024  Is_adjacent(RU,LUF)=false;
1025  Is_adjacent(RU,RDB)=false;
1026  Is_adjacent(RU,RDF)=false;
1027  Is_adjacent(RU,RUB)=true;
1028  Is_adjacent(RU,RUF)=true;
1029 
1030  Is_adjacent(RB,LDB)=false;
1031  Is_adjacent(RB,LDF)=false;
1032  Is_adjacent(RB,LUB)=false;
1033  Is_adjacent(RB,LUF)=false;
1034  Is_adjacent(RB,RDB)=true;
1035  Is_adjacent(RB,RDF)=false;
1036  Is_adjacent(RB,RUB)=true;
1037  Is_adjacent(RB,RUF)=false;
1038 
1039  Is_adjacent(RF,LDB)=false;
1040  Is_adjacent(RF,LDF)=false;
1041  Is_adjacent(RF,LUB)=false;
1042  Is_adjacent(RF,LUF)=false;
1043  Is_adjacent(RF,RDB)=false;
1044  Is_adjacent(RF,RDF)=true;
1045  Is_adjacent(RF,RUB)=false;
1046  Is_adjacent(RF,RUF)=true;
1047 
1048 
1049  Is_adjacent(DB,LDB)=true;
1050  Is_adjacent(DB,LDF)=false;
1051  Is_adjacent(DB,LUB)=false;
1052  Is_adjacent(DB,LUF)=false;
1053  Is_adjacent(DB,RDB)=true;
1054  Is_adjacent(DB,RDF)=false;
1055  Is_adjacent(DB,RUB)=false;
1056  Is_adjacent(DB,RUF)=false;
1057 
1058 
1059  Is_adjacent(DF,LDB)=false;
1060  Is_adjacent(DF,LDF)=true;
1061  Is_adjacent(DF,LUB)=false;
1062  Is_adjacent(DF,LUF)=false;
1063  Is_adjacent(DF,RDB)=false;
1064  Is_adjacent(DF,RDF)=true;
1065  Is_adjacent(DF,RUB)=false;
1066  Is_adjacent(DF,RUF)=false;
1067 
1068  Is_adjacent(UB,LDB)=false;
1069  Is_adjacent(UB,LDF)=false;
1070  Is_adjacent(UB,LUB)=true;
1071  Is_adjacent(UB,LUF)=false;
1072  Is_adjacent(UB,RDB)=false;
1073  Is_adjacent(UB,RDF)=false;
1074  Is_adjacent(UB,RUB)=true;
1075  Is_adjacent(UB,RUF)=false;
1076 
1077  Is_adjacent(UF,LDB)=false;
1078  Is_adjacent(UF,LDF)=false;
1079  Is_adjacent(UF,LUB)=false;
1080  Is_adjacent(UF,LUF)=true;
1081  Is_adjacent(UF,RDB)=false;
1082  Is_adjacent(UF,RDF)=false;
1083  Is_adjacent(UF,RUB)=false;
1084  Is_adjacent(UF,RUF)=true;
1085 
1086 
1087 
1088  // Common face of various octants from Samet's book
1089  Common_face.resize(27,27);
1093  Common_face(LDB,LUF)=L;
1095  Common_face(LDB,RDF)=D;
1096  Common_face(LDB,RUB)=B;
1098 
1101  Common_face(LDF,LUB)=L;
1103  Common_face(LDF,RDB)=D;
1106  Common_face(LDF,RUF)=F;
1107 
1109  Common_face(LUB,LDF)=L;
1112  Common_face(LUB,RDB)=B;
1115  Common_face(LUB,RUF)=U;
1116 
1117  Common_face(LUF,LDB)=L;
1122  Common_face(LUF,RDF)=F;
1123  Common_face(LUF,RUB)=U;
1125 
1127  Common_face(RDB,LDF)=D;
1128  Common_face(RDB,LUB)=B;
1133  Common_face(RDB,RUF)=R;
1134 
1135  Common_face(RDF,LDB)=D;
1138  Common_face(RDF,LUF)=F;
1141  Common_face(RDF,RUB)=R;
1143 
1144  Common_face(RUB,LDB)=B;
1147  Common_face(RUB,LUF)=U;
1149  Common_face(RUB,RDF)=R;
1152 
1154  Common_face(RUF,LDF)=F;
1155  Common_face(RUF,LUB)=U;
1157  Common_face(RUF,RDB)=R;
1161 
1162 
1163 
1164  // Common face of various edges/octants from Samet's book
1167  Common_face(LD,LUB)=L;
1168  Common_face(LD,LUF)=L;
1169  Common_face(LD,RDB)=D;
1170  Common_face(LD,RDF)=D;
1173 
1174  Common_face(LU,LDB)=L;
1175  Common_face(LU,LDF)=L;
1180  Common_face(LU,RUB)=U;
1181  Common_face(LU,RUF)=U;
1182 
1184  Common_face(LB,LDF)=L;
1186  Common_face(LB,LUF)=L;
1187  Common_face(LB,RDB)=B;
1189  Common_face(LB,RUB)=B;
1191 
1192  Common_face(LF,LDB)=L;
1194  Common_face(LF,LUB)=L;
1197  Common_face(LF,RDF)=F;
1199  Common_face(LF,RUF)=F;
1200 
1201  Common_face(RD,LDB)=D;
1202  Common_face(RD,LDF)=D;
1207  Common_face(RD,RUB)=R;
1208  Common_face(RD,RUF)=R;
1209 
1212  Common_face(RU,LUB)=U;
1213  Common_face(RU,LUF)=U;
1214  Common_face(RU,RDB)=R;
1215  Common_face(RU,RDF)=R;
1218 
1219  Common_face(RB,LDB)=B;
1221  Common_face(RB,LUB)=B;
1224  Common_face(RB,RDF)=R;
1226  Common_face(RB,RUF)=R;
1227 
1229  Common_face(RF,LDF)=F;
1231  Common_face(RF,LUF)=F;
1232  Common_face(RF,RDB)=R;
1234  Common_face(RF,RUB)=R;
1236 
1237 
1239  Common_face(DB,LDF)=D;
1240  Common_face(DB,LUB)=B;
1243  Common_face(DB,RDF)=D;
1244  Common_face(DB,RUB)=B;
1246 
1247  Common_face(DF,LDB)=D;
1250  Common_face(DF,LUF)=F;
1251  Common_face(DF,RDB)=D;
1254  Common_face(DF,RUF)=F;
1255 
1256  Common_face(UB,LDB)=B;
1259  Common_face(UB,LUF)=U;
1260  Common_face(UB,RDB)=B;
1263  Common_face(UB,RUF)=U;
1264 
1266  Common_face(UF,LDF)=F;
1267  Common_face(UF,LUB)=U;
1270  Common_face(UF,RDF)=F;
1271  Common_face(UF,RUB)=U;
1273 
1274 
1275  // Tecplot colours for neighbours in various directions
1276  Colour.resize(27);
1277  Colour[R]="CYAN";
1278  Colour[L]="RED";
1279  Colour[U]="GREEN";
1280  Colour[D]="BLUE";
1281  Colour[F]="CUSTOM3";
1282  Colour[B]="PURPLE";
1283  Colour[OMEGA]="YELLOW";
1284 
1285  Colour[LB]="BLUE";
1286  Colour[RB]="BLUE";
1287  Colour[DB]="BLUE";
1288  Colour[UB]="BLUE";
1289 
1290  Colour[LD]="GREEN";
1291  Colour[RD]="GREEN";
1292  Colour[LU]="GREEN";
1293  Colour[RU]="GREEN";
1294 
1295 
1296  Colour[LF]="RED";
1297  Colour[RF]="RED";
1298  Colour[DF]="RED";
1299  Colour[UF]="RED";
1300 
1301  Colour[OMEGA]="YELLOW";
1302 
1303 
1304 
1305  // Reflection scheme:
1306  // Reflect(direction,octant): Get mirror of octant in direction
1307  // Table in Samet book as well
1308  Reflect.resize(27,27);
1309 
1310  Reflect(L,LDB)=RDB;
1311  Reflect(R,LDB)=RDB;
1312  Reflect(D,LDB)=LUB;
1313  Reflect(U,LDB)=LUB;
1314  Reflect(B,LDB)=LDF;
1315  Reflect(F,LDB)=LDF;
1316 
1317  Reflect(L,LDF)=RDF;
1318  Reflect(R,LDF)=RDF;
1319  Reflect(D,LDF)=LUF;
1320  Reflect(U,LDF)=LUF;
1321  Reflect(B,LDF)=LDB;
1322  Reflect(F,LDF)=LDB;
1323 
1324  Reflect(L,LUB)=RUB;
1325  Reflect(R,LUB)=RUB;
1326  Reflect(D,LUB)=LDB;
1327  Reflect(U,LUB)=LDB;
1328  Reflect(B,LUB)=LUF;
1329  Reflect(F,LUB)=LUF;
1330 
1331  Reflect(L,LUF)=RUF;
1332  Reflect(R,LUF)=RUF;
1333  Reflect(D,LUF)=LDF;
1334  Reflect(U,LUF)=LDF;
1335  Reflect(B,LUF)=LUB;
1336  Reflect(F,LUF)=LUB;
1337 
1338  Reflect(L,RDB)=LDB;
1339  Reflect(R,RDB)=LDB;
1340  Reflect(D,RDB)=RUB;
1341  Reflect(U,RDB)=RUB;
1342  Reflect(B,RDB)=RDF;
1343  Reflect(F,RDB)=RDF;
1344 
1345  Reflect(L,RDF)=LDF;
1346  Reflect(R,RDF)=LDF;
1347  Reflect(D,RDF)=RUF;
1348  Reflect(U,RDF)=RUF;
1349  Reflect(B,RDF)=RDB;
1350  Reflect(F,RDF)=RDB;
1351 
1352  Reflect(L,RUB)=LUB;
1353  Reflect(R,RUB)=LUB;
1354  Reflect(D,RUB)=RDB;
1355  Reflect(U,RUB)=RDB;
1356  Reflect(B,RUB)=RUF;
1357  Reflect(F,RUB)=RUF;
1358 
1359  Reflect(L,RUF)=LUF;
1360  Reflect(R,RUF)=LUF;
1361  Reflect(D,RUF)=RDF;
1362  Reflect(U,RUF)=RDF;
1363  Reflect(B,RUF)=RUB;
1364  Reflect(F,RUF)=RUB;
1365 
1366 
1367  // Reflection scheme:
1368  // Reflect(direction (edge) ,octant): Get mirror of edge in direction
1369  // Table in Samet book as well
1370 
1371  Reflect(LD,LDB)=RUB;
1372  Reflect(LU,LDB)=RUB;
1373  Reflect(RD,LDB)=RUB;
1374  Reflect(RU,LDB)=RUB;
1375 
1376  Reflect(LD,LDF)=RUF;
1377  Reflect(LU,LDF)=RUF;
1378  Reflect(RD,LDF)=RUF;
1379  Reflect(RU,LDF)=RUF;
1380 
1381  Reflect(LD,LUB)=RDB;
1382  Reflect(LU,LUB)=RDB;
1383  Reflect(RD,LUB)=RDB;
1384  Reflect(RU,LUB)=RDB;
1385 
1386  Reflect(LD,LUF)=RDF;
1387  Reflect(LU,LUF)=RDF;
1388  Reflect(RD,LUF)=RDF;
1389  Reflect(RU,LUF)=RDF;
1390 
1391 
1392  Reflect(LD,RDB)=LUB;
1393  Reflect(LU,RDB)=LUB;
1394  Reflect(RD,RDB)=LUB;
1395  Reflect(RU,RDB)=LUB;
1396 
1397  Reflect(LD,RDF)=LUF;
1398  Reflect(LU,RDF)=LUF;
1399  Reflect(RD,RDF)=LUF;
1400  Reflect(RU,RDF)=LUF;
1401 
1402  Reflect(LD,RUB)=LDB;
1403  Reflect(LU,RUB)=LDB;
1404  Reflect(RD,RUB)=LDB;
1405  Reflect(RU,RUB)=LDB;
1406 
1407  Reflect(LD,RUF)=LDF;
1408  Reflect(LU,RUF)=LDF;
1409  Reflect(RD,RUF)=LDF;
1410  Reflect(RU,RUF)=LDF;
1411 
1412 
1413 
1414  Reflect(LB,LDB)=RDF;
1415  Reflect(LF,LDB)=RDF;
1416  Reflect(RB,LDB)=RDF;
1417  Reflect(RF,LDB)=RDF;
1418 
1419  Reflect(LB,LDF)=RDB;
1420  Reflect(LF,LDF)=RDB;
1421  Reflect(RB,LDF)=RDB;
1422  Reflect(RF,LDF)=RDB;
1423 
1424  Reflect(LB,LUB)=RUF;
1425  Reflect(LF,LUB)=RUF;
1426  Reflect(RB,LUB)=RUF;
1427  Reflect(RF,LUB)=RUF;
1428 
1429  Reflect(LB,LUF)=RUB;
1430  Reflect(LF,LUF)=RUB;
1431  Reflect(RB,LUF)=RUB;
1432  Reflect(RF,LUF)=RUB;
1433 
1434  Reflect(LB,RDB)=LDF;
1435  Reflect(LF,RDB)=LDF;
1436  Reflect(RB,RDB)=LDF;
1437  Reflect(RF,RDB)=LDF;
1438 
1439  Reflect(LB,RDF)=LDB;
1440  Reflect(LF,RDF)=LDB;
1441  Reflect(RB,RDF)=LDB;
1442  Reflect(RF,RDF)=LDB;
1443 
1444  Reflect(LB,RUB)=LUF;
1445  Reflect(LF,RUB)=LUF;
1446  Reflect(RB,RUB)=LUF;
1447  Reflect(RF,RUB)=LUF;
1448 
1449  Reflect(LB,RUF)=LUB;
1450  Reflect(LF,RUF)=LUB;
1451  Reflect(RB,RUF)=LUB;
1452  Reflect(RF,RUF)=LUB;
1453 
1454 
1455  Reflect(DB,LDB)=LUF;
1456  Reflect(DF,LDB)=LUF;
1457  Reflect(UB,LDB)=LUF;
1458  Reflect(UF,LDB)=LUF;
1459 
1460  Reflect(DB,LDF)=LUB;
1461  Reflect(DF,LDF)=LUB;
1462  Reflect(UB,LDF)=LUB;
1463  Reflect(UF,LDF)=LUB;
1464 
1465  Reflect(DB,LUB)=LDF;
1466  Reflect(DF,LUB)=LDF;
1467  Reflect(UB,LUB)=LDF;
1468  Reflect(UF,LUB)=LDF;
1469 
1470  Reflect(DB,LUF)=LDB;
1471  Reflect(DF,LUF)=LDB;
1472  Reflect(UB,LUF)=LDB;
1473  Reflect(UF,LUF)=LDB;
1474 
1475  Reflect(DB,RDB)=RUF;
1476  Reflect(DF,RDB)=RUF;
1477  Reflect(UB,RDB)=RUF;
1478  Reflect(UF,RDB)=RUF;
1479 
1480  Reflect(DB,RDF)=RUB;
1481  Reflect(DF,RDF)=RUB;
1482  Reflect(UB,RDF)=RUB;
1483  Reflect(UF,RDF)=RUB;
1484 
1485  Reflect(DB,RUB)=RDF;
1486  Reflect(DF,RUB)=RDF;
1487  Reflect(UB,RUB)=RDF;
1488  Reflect(UF,RUB)=RDF;
1489 
1490  Reflect(DB,RUF)=RDB;
1491  Reflect(DF,RUF)=RDB;
1492  Reflect(UB,RUF)=RDB;
1493  Reflect(UF,RUF)=RDB;
1494 
1495 
1496 
1497  // S_base(i,direction) etc.: Initial value/increment for coordinate s[i] on
1498  // the face indicated by direction (L/R/U/D/F/B)
1499  S_base.resize(3,27);
1500  S_steplo.resize(3,27);
1501  S_stephi.resize(3,27);
1502 
1503  S_base(0,L)=-1.0;
1504  S_base(1,L)=-1.0;
1505  S_base(2,L)=-1.0;
1506  S_steplo(0,L)=0.0;
1507  S_steplo(1,L)=2.0;
1508  S_steplo(2,L)=0.0;
1509  S_stephi(0,L)=0.0;
1510  S_stephi(1,L)=0.0;
1511  S_stephi(2,L)=2.0;
1512 
1513  S_base(0,R)=1.0;
1514  S_base(1,R)=-1.0;
1515  S_base(2,R)=-1.0;
1516  S_steplo(0,R)=0.0;
1517  S_steplo(1,R)=2.0;
1518  S_steplo(2,R)=0.0;
1519  S_stephi(0,R)=0.0;
1520  S_stephi(1,R)=0.0;
1521  S_stephi(2,R)=2.0;
1522 
1523  S_base(0,U)=-1.0;
1524  S_base(1,U)=1.0;
1525  S_base(2,U)=-1.0;
1526  S_steplo(0,U)=2.0;
1527  S_steplo(1,U)=0.0;
1528  S_steplo(2,U)=0.0;
1529  S_stephi(0,U)=0.0;
1530  S_stephi(1,U)=0.0;
1531  S_stephi(2,U)=2.0;
1532 
1533  S_base(0,D)=-1.0;
1534  S_base(1,D)=-1.0;
1535  S_base(2,D)=-1.0;
1536  S_steplo(0,D)=2.0;
1537  S_steplo(1,D)=0.0;
1538  S_steplo(2,D)=0.0;
1539  S_stephi(0,D)=0.0;
1540  S_stephi(1,D)=0.0;
1541  S_stephi(2,D)=2.0;
1542 
1543  S_base(0,F)=-1.0;
1544  S_base(1,F)=-1.0;
1545  S_base(2,F)=1.0;
1546  S_steplo(0,F)=2.0;
1547  S_steplo(1,F)=0.0;
1548  S_steplo(2,F)=0.0;
1549  S_stephi(0,F)=0.0;
1550  S_stephi(1,F)=2.0;
1551  S_stephi(2,F)=0.0;
1552 
1553  S_base(0,B)=-1.0;
1554  S_base(1,B)=-1.0;
1555  S_base(2,B)=-1.0;
1556  S_steplo(0,B)=2.0;
1557  S_steplo(1,B)=0.0;
1558  S_steplo(2,B)=0.0;
1559  S_stephi(0,B)=0.0;
1560  S_stephi(1,B)=2.0;
1561  S_stephi(2,B)=0.0;
1562 
1563 
1564  // Relative to the left/down/back vertex in any (father) octree, the
1565  // corresponding vertex in the son specified by \c son_octant has an offset.
1566  // If we project the son_octant's left/down/back vertex onto the
1567  // father's face \c face, it is located at the in-face coordinate
1568  // \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
1569  // \c S_steplo for an explanation of the subscripts \c _hi and \c _lo.]
1570  S_directlo.resize(27,27);
1571  S_directhi.resize(27,27);
1572 
1573  S_directlo(L,LDB)=0.0;
1574  S_directlo(R,LDB)=0.0;
1575  S_directlo(U,LDB)=0.0;
1576  S_directlo(D,LDB)=0.0;
1577  S_directlo(F,LDB)=0.0;
1578  S_directlo(B,LDB)=0.0;
1579 
1580  S_directlo(L,LDF)=0.0;
1581  S_directlo(R,LDF)=0.0;
1582  S_directlo(U,LDF)=0.0;
1583  S_directlo(D,LDF)=0.0;
1584  S_directlo(F,LDF)=0.0;
1585  S_directlo(B,LDF)=0.0;
1586 
1587  S_directlo(L,LUB)=1.0;
1588  S_directlo(R,LUB)=1.0;
1589  S_directlo(U,LUB)=0.0;
1590  S_directlo(D,LUB)=0.0;
1591  S_directlo(F,LUB)=0.0;
1592  S_directlo(B,LUB)=0.0;
1593 
1594  S_directlo(L,LUF)=1.0;
1595  S_directlo(R,LUF)=1.0;
1596  S_directlo(U,LUF)=0.0;
1597  S_directlo(D,LUF)=0.0;
1598  S_directlo(F,LUF)=0.0;
1599  S_directlo(B,LUF)=0.0;
1600 
1601  S_directlo(L,RDB)=0.0;
1602  S_directlo(R,RDB)=0.0;
1603  S_directlo(U,RDB)=1.0;
1604  S_directlo(D,RDB)=1.0;
1605  S_directlo(F,RDB)=1.0;
1606  S_directlo(B,RDB)=1.0;
1607 
1608  S_directlo(L,RDF)=0.0;
1609  S_directlo(R,RDF)=0.0;
1610  S_directlo(U,RDF)=1.0;
1611  S_directlo(D,RDF)=1.0;
1612  S_directlo(F,RDF)=1.0;
1613  S_directlo(B,RDF)=1.0;
1614 
1615  S_directlo(L,RUB)=1.0;
1616  S_directlo(R,RUB)=1.0;
1617  S_directlo(U,RUB)=1.0;
1618  S_directlo(D,RUB)=1.0;
1619  S_directlo(F,RUB)=1.0;
1620  S_directlo(B,RUB)=1.0;
1621 
1622  S_directlo(L,RUF)=1.0;
1623  S_directlo(R,RUF)=1.0;
1624  S_directlo(U,RUF)=1.0;
1625  S_directlo(D,RUF)=1.0;
1626  S_directlo(F,RUF)=1.0;
1627  S_directlo(B,RUF)=1.0;
1628 
1629 
1630 
1631 
1632  S_directhi(L,LDB)=0.0;
1633  S_directhi(R,LDB)=0.0;
1634  S_directhi(U,LDB)=0.0;
1635  S_directhi(D,LDB)=0.0;
1636  S_directhi(F,LDB)=0.0;
1637  S_directhi(B,LDB)=0.0;
1638 
1639  S_directhi(L,LDF)=1.0;
1640  S_directhi(R,LDF)=1.0;
1641  S_directhi(U,LDF)=1.0;
1642  S_directhi(D,LDF)=1.0;
1643  S_directhi(F,LDF)=0.0;
1644  S_directhi(B,LDF)=0.0;
1645 
1646  S_directhi(L,LUB)=0.0;
1647  S_directhi(R,LUB)=0.0;
1648  S_directhi(U,LUB)=0.0;
1649  S_directhi(D,LUB)=0.0;
1650  S_directhi(F,LUB)=1.0;
1651  S_directhi(B,LUB)=1.0;
1652 
1653  S_directhi(L,LUF)=1.0;
1654  S_directhi(R,LUF)=1.0;
1655  S_directhi(U,LUF)=1.0;
1656  S_directhi(D,LUF)=1.0;
1657  S_directhi(F,LUF)=1.0;
1658  S_directhi(B,LUF)=1.0;
1659 
1660  S_directhi(L,RDB)=0.0;
1661  S_directhi(R,RDB)=0.0;
1662  S_directhi(U,RDB)=0.0;
1663  S_directhi(D,RDB)=0.0;
1664  S_directhi(F,RDB)=0.0;
1665  S_directhi(B,RDB)=0.0;
1666 
1667  S_directhi(L,RDF)=1.0;
1668  S_directhi(R,RDF)=1.0;
1669  S_directhi(U,RDF)=1.0;
1670  S_directhi(D,RDF)=1.0;
1671  S_directhi(F,RDF)=0.0;
1672  S_directhi(B,RDF)=0.0;
1673 
1674  S_directhi(L,RUB)=0.0;
1675  S_directhi(R,RUB)=0.0;
1676  S_directhi(U,RUB)=0.0;
1677  S_directhi(D,RUB)=0.0;
1678  S_directhi(F,RUB)=1.0;
1679  S_directhi(B,RUB)=1.0;
1680 
1681  S_directhi(L,RUF)=1.0;
1682  S_directhi(R,RUF)=1.0;
1683  S_directhi(U,RUF)=1.0;
1684  S_directhi(D,RUF)=1.0;
1685  S_directhi(F,RUF)=1.0;
1686  S_directhi(B,RUF)=1.0;
1687 
1688 
1689 
1690 
1691 
1692 
1693 
1694  // S_base_edge(i,direction): Initial value/increment for coordinate s[i] on
1695  // the edge indicated by direction (LB,RB,...)
1696  S_base_edge.resize(3,27);
1697  S_step_edge.resize(3,27);
1698 
1699  S_base_edge(0,LB)=-1.0;
1700  S_base_edge(1,LB)=-1.0;
1701  S_base_edge(2,LB)=-1.0;
1702  S_step_edge(0,LB)=0.0;
1703  S_step_edge(1,LB)=2.0;
1704  S_step_edge(2,LB)=0.0;
1705 
1706  S_base_edge(0,RB)= 1.0;
1707  S_base_edge(1,RB)=-1.0;
1708  S_base_edge(2,RB)=-1.0;
1709  S_step_edge(0,RB)=0.0;
1710  S_step_edge(1,RB)=2.0;
1711  S_step_edge(2,RB)=0.0;
1712 
1713  S_base_edge(0,DB)=-1.0;
1714  S_base_edge(1,DB)=-1.0;
1715  S_base_edge(2,DB)=-1.0;
1716  S_step_edge(0,DB)=2.0;
1717  S_step_edge(1,DB)=0.0;
1718  S_step_edge(2,DB)=0.0;
1719 
1720  S_base_edge(0,UB)=-1.0;
1721  S_base_edge(1,UB)= 1.0;
1722  S_base_edge(2,UB)=-1.0;
1723  S_step_edge(0,UB)=2.0;
1724  S_step_edge(1,UB)=0.0;
1725  S_step_edge(2,UB)=0.0;
1726 
1727  S_base_edge(0,LD)=-1.0;
1728  S_base_edge(1,LD)=-1.0;
1729  S_base_edge(2,LD)=-1.0;
1730  S_step_edge(0,LD)=0.0;
1731  S_step_edge(1,LD)=0.0;
1732  S_step_edge(2,LD)=2.0;
1733 
1734  S_base_edge(0,RD)= 1.0;
1735  S_base_edge(1,RD)=-1.0;
1736  S_base_edge(2,RD)=-1.0;
1737  S_step_edge(0,RD)=0.0;
1738  S_step_edge(1,RD)=0.0;
1739  S_step_edge(2,RD)=2.0;
1740 
1741  S_base_edge(0,LU)=-1.0;
1742  S_base_edge(1,LU)= 1.0;
1743  S_base_edge(2,LU)=-1.0;
1744  S_step_edge(0,LU)=0.0;
1745  S_step_edge(1,LU)=0.0;
1746  S_step_edge(2,LU)=2.0;
1747 
1748  S_base_edge(0,RU)= 1.0;
1749  S_base_edge(1,RU)= 1.0;
1750  S_base_edge(2,RU)=-1.0;
1751  S_step_edge(0,RU)=0.0;
1752  S_step_edge(1,RU)=0.0;
1753  S_step_edge(2,RU)=2.0;
1754 
1755 
1756 
1757  S_base_edge(0,LF)=-1.0;
1758  S_base_edge(1,LF)=-1.0;
1759  S_base_edge(2,LF)= 1.0;
1760  S_step_edge(0,LF)=0.0;
1761  S_step_edge(1,LF)=2.0;
1762  S_step_edge(2,LF)=0.0;
1763 
1764  S_base_edge(0,RF)= 1.0;
1765  S_base_edge(1,RF)=-1.0;
1766  S_base_edge(2,RF)= 1.0;
1767  S_step_edge(0,RF)=0.0;
1768  S_step_edge(1,RF)=2.0;
1769  S_step_edge(2,RF)=0.0;
1770 
1771  S_base_edge(0,DF)=-1.0;
1772  S_base_edge(1,DF)=-1.0;
1773  S_base_edge(2,DF)= 1.0;
1774  S_step_edge(0,DF)=2.0;
1775  S_step_edge(1,DF)=0.0;
1776  S_step_edge(2,DF)=0.0;
1777 
1778  S_base_edge(0,UF)=-1.0;
1779  S_base_edge(1,UF)= 1.0;
1780  S_base_edge(2,UF)= 1.0;
1781  S_step_edge(0,UF)=2.0;
1782  S_step_edge(1,UF)=0.0;
1783  S_step_edge(2,UF)=0.0;
1784 
1785 
1786 
1787 
1788  // Relative to the left/down/back vertex in any (father) octree, the
1789  // corresponding vertex in the son specified by \c son_octant has an offset.
1790  // If we project the son_octant's left/down/back vertex onto the
1791  // father's edge \c edge, it is located at the in-face coordinate
1792  // \c s_lo = h/2 \c S_direct_edge(edge,son_octant).
1793  S_direct_edge.resize(27,27);
1794  S_direct_edge(LB,LDB)=0.0;
1795  S_direct_edge(RB,LDB)=0.0;
1796  S_direct_edge(DB,LDB)=0.0;
1797  S_direct_edge(UB,LDB)=0.0;
1798  S_direct_edge(LD,LDB)=0.0;
1799  S_direct_edge(RD,LDB)=0.0;
1800  S_direct_edge(LU,LDB)=0.0;
1801  S_direct_edge(RU,LDB)=0.0;
1802  S_direct_edge(LF,LDB)=0.0;
1803  S_direct_edge(RF,LDB)=0.0;
1804  S_direct_edge(DF,LDB)=0.0;
1805  S_direct_edge(UF,LDB)=0.0;
1806 
1807  S_direct_edge(LB,RDB)=0.0;
1808  S_direct_edge(RB,RDB)=0.0;
1809  S_direct_edge(DB,RDB)=1.0;
1810  S_direct_edge(UB,RDB)=1.0;
1811  S_direct_edge(LD,RDB)=0.0;
1812  S_direct_edge(RD,RDB)=0.0;
1813  S_direct_edge(LU,RDB)=0.0;
1814  S_direct_edge(RU,RDB)=0.0;
1815  S_direct_edge(LF,RDB)=0.0;
1816  S_direct_edge(RF,RDB)=0.0;
1817  S_direct_edge(DF,RDB)=1.0;
1818  S_direct_edge(UF,RDB)=1.0;
1819 
1820  S_direct_edge(LB,LUB)=1.0;
1821  S_direct_edge(RB,LUB)=1.0;
1822  S_direct_edge(DB,LUB)=0.0;
1823  S_direct_edge(UB,LUB)=0.0;
1824  S_direct_edge(LD,LUB)=0.0;
1825  S_direct_edge(RD,LUB)=0.0;
1826  S_direct_edge(LU,LUB)=0.0;
1827  S_direct_edge(RU,LUB)=0.0;
1828  S_direct_edge(LF,LUB)=1.0;
1829  S_direct_edge(RF,LUB)=1.0;
1830  S_direct_edge(DF,LUB)=0.0;
1831  S_direct_edge(UF,LUB)=0.0;
1832 
1833  S_direct_edge(LB,RUB)=1.0;
1834  S_direct_edge(RB,RUB)=1.0;
1835  S_direct_edge(DB,RUB)=1.0;
1836  S_direct_edge(UB,RUB)=1.0;
1837  S_direct_edge(LD,RUB)=0.0;
1838  S_direct_edge(RD,RUB)=0.0;
1839  S_direct_edge(LU,RUB)=0.0;
1840  S_direct_edge(RU,RUB)=0.0;
1841  S_direct_edge(LF,RUB)=1.0;
1842  S_direct_edge(RF,RUB)=1.0;
1843  S_direct_edge(DF,RUB)=1.0;
1844  S_direct_edge(UF,RUB)=1.0;
1845 
1846 
1847  S_direct_edge(LB,LDF)=0.0;
1848  S_direct_edge(RB,LDF)=0.0;
1849  S_direct_edge(DB,LDF)=0.0;
1850  S_direct_edge(UB,LDF)=0.0;
1851  S_direct_edge(LD,LDF)=1.0;
1852  S_direct_edge(RD,LDF)=1.0;
1853  S_direct_edge(LU,LDF)=1.0;
1854  S_direct_edge(RU,LDF)=1.0;
1855  S_direct_edge(LF,LDF)=0.0;
1856  S_direct_edge(RF,LDF)=0.0;
1857  S_direct_edge(DF,LDF)=0.0;
1858  S_direct_edge(UF,LDF)=0.0;
1859 
1860  S_direct_edge(LB,RDF)=0.0;
1861  S_direct_edge(RB,RDF)=0.0;
1862  S_direct_edge(DB,RDF)=1.0;
1863  S_direct_edge(UB,RDF)=1.0;
1864  S_direct_edge(LD,RDF)=1.0;
1865  S_direct_edge(RD,RDF)=1.0;
1866  S_direct_edge(LU,RDF)=1.0;
1867  S_direct_edge(RU,RDF)=1.0;
1868  S_direct_edge(LF,RDF)=0.0;
1869  S_direct_edge(RF,RDF)=0.0;
1870  S_direct_edge(DF,RDF)=1.0;
1871  S_direct_edge(UF,RDF)=1.0;
1872 
1873  S_direct_edge(LB,LUF)=1.0;
1874  S_direct_edge(RB,LUF)=1.0;
1875  S_direct_edge(DB,LUF)=0.0;
1876  S_direct_edge(UB,LUF)=0.0;
1877  S_direct_edge(LD,LUF)=1.0;
1878  S_direct_edge(RD,LUF)=1.0;
1879  S_direct_edge(LU,LUF)=1.0;
1880  S_direct_edge(RU,LUF)=1.0;
1881  S_direct_edge(LF,LUF)=1.0;
1882  S_direct_edge(RF,LUF)=1.0;
1883  S_direct_edge(DF,LUF)=0.0;
1884  S_direct_edge(UF,LUF)=0.0;
1885 
1886  S_direct_edge(LB,RUF)=1.0;
1887  S_direct_edge(RB,RUF)=1.0;
1888  S_direct_edge(DB,RUF)=1.0;
1889  S_direct_edge(UB,RUF)=1.0;
1890  S_direct_edge(LD,RUF)=1.0;
1891  S_direct_edge(RD,RUF)=1.0;
1892  S_direct_edge(LU,RUF)=1.0;
1893  S_direct_edge(RU,RUF)=1.0;
1894  S_direct_edge(LF,RUF)=1.0;
1895  S_direct_edge(RF,RUF)=1.0;
1896  S_direct_edge(DF,RUF)=1.0;
1897  S_direct_edge(UF,RUF)=1.0;
1898 
1899 
1900 
1901 //Public static data:
1902 //-------------------
1903 
1904  // Translate (enumerated) directions into strings
1905  Direct_string.resize(27);
1906  Direct_string[LDB]="LDB";
1907  Direct_string[LDF]="LDF";
1908  Direct_string[LUB]="LUB";
1909  Direct_string[LUF]="LUF";
1910  Direct_string[RDB]="RDB";
1911  Direct_string[RDF]="RDF";
1912  Direct_string[RUB]="RUB";
1913  Direct_string[RUF]="RUF";
1914 
1915 
1916  Direct_string[L]="L";
1917  Direct_string[R]="R";
1918  Direct_string[U]="U";
1919  Direct_string[D]="D";
1920  Direct_string[F]="F";
1921  Direct_string[B]="B";
1922 
1923  Direct_string[LU]="LU";
1924  Direct_string[LD]="LD";
1925  Direct_string[LF]="LF";
1926  Direct_string[LB]="LB";
1927  Direct_string[RU]="RU";
1928  Direct_string[RD]="RD";
1929  Direct_string[RF]="RF";
1930  Direct_string[RB]="RB";
1931  Direct_string[UF]="UF";
1932  Direct_string[UB]="UB";
1933  Direct_string[DF]="DF";
1934  Direct_string[DB]="DB";
1935 
1936  Direct_string[OMEGA]="OMEGA";
1937 
1938 
1939  // Get opposite face, e.g. Reflect_face(L)=R
1940  Reflect_face.resize(27);
1941  Reflect_face[L]=R;
1942  Reflect_face[R]=L;
1943  Reflect_face[U]=D;
1944  Reflect_face[D]=U;
1945  Reflect_face[B]=F;
1946  Reflect_face[F]=B;
1947 
1948  // Get opposite edge, e.g. Reflect_edge(DB)=UF
1949  Reflect_edge.resize(27);
1950  Reflect_edge[LB]=RF;
1951  Reflect_edge[RB]=LF;
1952  Reflect_edge[DB]=UF;
1953  Reflect_edge[UB]=DF;
1954 
1955  Reflect_edge[LD]=RU;
1956  Reflect_edge[RD]=LU;
1957  Reflect_edge[LU]=RD;
1958  Reflect_edge[RU]=LD;
1959 
1960  Reflect_edge[LF]=RB;
1961  Reflect_edge[RF]=LB;
1962  Reflect_edge[DF]=UB;
1963  Reflect_edge[UF]=DB;
1964 
1965  // Get opposite vertex, e.g. Reflect_vertex(LDB)=RUF
1966  Reflect_vertex.resize(27);
1975 
1976 
1977 
1978  // Vertices at ends of edges
1979  Vertex_at_end_of_edge.resize(27);
1980 
1981  Vertex_at_end_of_edge[DB].resize(2);
1982  Vertex_at_end_of_edge[DB][0]=LDB; // Pattern: both other indices
1984 
1985  Vertex_at_end_of_edge[UB].resize(2);
1986  Vertex_at_end_of_edge[UB][0]=LUB; // Pattern: both other indices
1988 
1989 
1990  Vertex_at_end_of_edge[LB].resize(2);
1991  Vertex_at_end_of_edge[LB][0]=LUB; // Pattern: both other indices
1993 
1994  Vertex_at_end_of_edge[RB].resize(2);
1995  Vertex_at_end_of_edge[RB][0]=RUB; // Pattern: both other indices
1997 
1998 
1999 
2000  Vertex_at_end_of_edge[LD].resize(2);
2001  Vertex_at_end_of_edge[LD][0]=LDF; // Pattern: both other indices
2003 
2004  Vertex_at_end_of_edge[RD].resize(2);
2005  Vertex_at_end_of_edge[RD][0]=RDF; // Pattern: both other indices
2007 
2008  Vertex_at_end_of_edge[LU].resize(2);
2009  Vertex_at_end_of_edge[LU][0]=LUF; // Pattern: both other indices
2011 
2012  Vertex_at_end_of_edge[RU].resize(2);
2013  Vertex_at_end_of_edge[RU][0]=RUF; // Pattern: both other indices
2015 
2016 
2017 
2018  Vertex_at_end_of_edge[DF].resize(2);
2019  Vertex_at_end_of_edge[DF][0]=LDF; // Pattern: both other indices
2021 
2022  Vertex_at_end_of_edge[UF].resize(2);
2023  Vertex_at_end_of_edge[UF][0]=LUF; // Pattern: both other indices
2025 
2026  Vertex_at_end_of_edge[LF].resize(2);
2027  Vertex_at_end_of_edge[LF][0]=LUF; // Pattern: both other indices
2029 
2030  Vertex_at_end_of_edge[RF].resize(2);
2031  Vertex_at_end_of_edge[RF][0]=RUF; // Pattern: both other indices
2033 
2034 
2035 
2036 
2037  // Initialisation of the values of Vector_to_direction
2038  Vector<int> vect(3);
2039  int elem;
2040 
2041  for(int i=-1; i<2; i++)
2042  {
2043  for(int j=-1; j<2; j++)
2044  {
2045  for(int k=-1; k<2; k++)
2046  {
2047  vect[0]=i; vect[1]=j; vect[2]=k;
2048  int num_elem=0;
2049 
2050  // To put a number on the vector (i,j,k), we assume that that
2051  // the vector (i+1,j+1,k+1) represents the decomposition
2052  // of the number of the corresponding direction in base 3.
2053  num_elem=(i+1)*9+(j+1)*3+(k+1);
2054 
2055  //for each number we have the corresponding element
2056  switch(num_elem)
2057  {
2058  case 6 :elem=LUB; break;
2059  case 24 :elem=RUB; break;
2060  case 26 :elem=RUF; break;
2061  case 8 :elem=LUF; break;
2062  case 0 :elem=LDB; break;
2063  case 18 :elem=RDB; break;
2064  case 20 :elem=RDF; break;
2065  case 2 :elem=LDF; break;
2066  case 25 :elem=RU; break;
2067  case 23 :elem=RF; break;
2068  case 19 :elem=RD; break;
2069  case 21 :elem=RB; break;
2070  case 7 :elem=LU; break;
2071  case 5 :elem=LF; break;
2072  case 1 :elem=LD; break;
2073  case 3 :elem=LB; break;
2074  case 17 :elem=UF; break;
2075  case 15 :elem=UB; break;
2076  case 11 :elem=DF; break;
2077  case 9 :elem=DB; break;
2078  case 16 :elem=U; break;
2079  case 10 :elem=D; break;
2080  case 22 :elem=R; break;
2081  case 4 :elem=L; break;
2082  case 14 :elem=F; break;
2083  case 12 :elem=B; break;
2084  case 13 :elem=OMEGA; break;
2085  default :
2086  elem=OMEGA;
2087  oomph_info
2088  << "there might be a problem with Vector_to_direction"
2089  <<std::endl;
2090  break;
2091  }
2092  Vector_to_direction[vect]=elem;
2093  }
2094  }
2095  }
2096 
2097 
2098 
2099  // Initialisation of Direction_to_vector:
2100  // Translate Octant, face, edge into direction vector using the
2101  // value of Direct_string; Direction_to_vector[U]={0,1,0};
2102  Direction_to_vector.resize(27);
2103  for(int i=LDB;i<=F;i++)
2104  {
2105  Direction_to_vector[i].resize(3);
2106  // Initialisation to 0;
2107  for(int j=0;j<3;j++)
2108  {
2109  Direction_to_vector[i][j]=0;
2110  }
2111 
2112  // Use +1 or -1 to indicate the relevant components of
2113  // the vector Direction
2114  for (unsigned j = 0; j < 3; j++)
2115  {
2116  if (Direct_string[i].length() > j)
2117  {
2118  switch (Direct_string[i][j])
2119  {
2120  case 'R':
2121  Direction_to_vector[i][0] = 1;
2122  break;
2123  case 'L':
2124  Direction_to_vector[i][0] = -1;
2125  break;
2126  case 'U':
2127  Direction_to_vector[i][1] = 1;
2128  break;
2129  case 'D':
2130  Direction_to_vector[i][1] = -1;
2131  break;
2132  case 'F':
2133  Direction_to_vector[i][2] = 1;
2134  break;
2135  case 'B':
2136  Direction_to_vector[i][2] = -1;
2137  break;
2138  default:
2139  oomph_info << "Direction Error !!" << std::endl;
2140  }
2141  }
2142  }
2143  }
2144 
2145 
2146  // Setup map that works out required rotations based on
2147  //-----------------------------------------------------
2148  // adjacent edge vertices
2149  //-----------------------
2150 
2151  int new_up,new_right;
2152  int new_vertex;
2153 
2154 
2155  // Map that stores the set of rotations (as a pairs consisting of
2156  // the up_equivalent and the right_equivalent) that move the vertex specified
2157  // by the first entry in key pair to the position of the second one:
2158  std::map<std::pair<int,int>, std::set<std::pair<int,int> > >
2159  required_rotation;
2160 
2161  // Loop over all vertices
2162  for (int vertex=LDB;vertex<=RUF;vertex++)
2163  {
2164 
2165  new_up=U;
2166  new_right=R;
2167  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2168  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2169  std::make_pair(new_up,new_right));
2170 
2171  new_up=U;
2172  new_right=B;
2173  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2174  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2175  std::make_pair(new_up,new_right));
2176 
2177  new_up=U;
2178  new_right=L;
2179  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2180  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2181  std::make_pair(new_up,new_right));
2182 
2183  new_up=U;
2184  new_right=F;
2185  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2186  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2187  std::make_pair(new_up,new_right));
2188 
2189 
2190 
2191 
2192 
2193 
2194  new_up=D;
2195  new_right=R;
2196  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2197  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2198  std::make_pair(new_up,new_right));
2199 
2200  new_up=D;
2201  new_right=B;
2202  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2203  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2204  std::make_pair(new_up,new_right));
2205 
2206  new_up=D;
2207  new_right=L;
2208  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2209  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2210  std::make_pair(new_up,new_right));
2211 
2212  new_up=D;
2213  new_right=F;
2214  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2215  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2216  std::make_pair(new_up,new_right));
2217 
2218 
2219 
2220 
2221  new_up=R;
2222  new_right=D;
2223  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2224  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2225  std::make_pair(new_up,new_right));
2226 
2227  new_up=R;
2228  new_right=B;
2229  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2230  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2231  std::make_pair(new_up,new_right));
2232 
2233  new_up=R;
2234  new_right=U;
2235  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2236  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2237  std::make_pair(new_up,new_right));
2238 
2239  new_up=R;
2240  new_right=F;
2241  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2242  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2243  std::make_pair(new_up,new_right));
2244 
2245 
2246 
2247 
2248 
2249  new_up=L;
2250  new_right=D;
2251  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2252  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2253  std::make_pair(new_up,new_right));
2254 
2255  new_up=L;
2256  new_right=B;
2257  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2258  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2259  std::make_pair(new_up,new_right));
2260 
2261  new_up=L;
2262  new_right=U;
2263  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2264  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2265  std::make_pair(new_up,new_right));
2266 
2267  new_up=L;
2268  new_right=F;
2269  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2270  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2271  std::make_pair(new_up,new_right));
2272 
2273 
2274 
2275  new_up=F;
2276  new_right=R;
2277  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2278  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2279  std::make_pair(new_up,new_right));
2280 
2281  new_up=F;
2282  new_right=L;
2283  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2284  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2285  std::make_pair(new_up,new_right));
2286 
2287  new_up=F;
2288  new_right=U;
2289  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2290  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2291  std::make_pair(new_up,new_right));
2292 
2293  new_up=F;
2294  new_right=D;
2295  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2296  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2297  std::make_pair(new_up,new_right));
2298 
2299 
2300 
2301 
2302 
2303 
2304  new_up=B;
2305  new_right=R;
2306  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2307  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2308  std::make_pair(new_up,new_right));
2309 
2310  new_up=B;
2311  new_right=L;
2312  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2313  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2314  std::make_pair(new_up,new_right));
2315 
2316  new_up=B;
2317  new_right=U;
2318  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2319  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2320  std::make_pair(new_up,new_right));
2321 
2322  new_up=B;
2323  new_right=D;
2324  new_vertex=OcTree::rotate(new_up,new_right,vertex);
2325  required_rotation[std::make_pair(vertex,new_vertex)].insert(
2326  std::make_pair(new_up,new_right));
2327 
2328 
2329  }
2330 
2331 
2332  // Each vertex is part of three edges. This container stores the
2333  // vertices in each of the three edge neighbours that are
2334  // adjacent to this node if there's no relative rotation between
2335  // the elements.
2336  std::map<int,Vector<int> > vertex_in_edge_neighbour;
2337 
2338 
2339  // Each vertex is part of three edges. This container stores the
2340  // vertices at the other end of the edge
2341  std::map<int,Vector<int> > other_vertex_on_edge;
2342 
2343  // Each vertex is part of three edges. This container stores the
2344  // vertices in the adjacent element at the other end of the edge
2345  // assuming there are no rotations between the elements.
2346  std::map<int,Vector<int> > other_vertex_in_edge_neighbour;
2347 
2348 
2349 
2350  vertex_in_edge_neighbour[LDB].resize(3);
2351  vertex_in_edge_neighbour[LDB][0]=LUF; // Pattern: exactly one letter matches
2352  vertex_in_edge_neighbour[LDB][1]=RDF;
2353  vertex_in_edge_neighbour[LDB][2]=RUB;
2354 
2355  other_vertex_on_edge[LDB].resize(3);
2356  other_vertex_on_edge[LDB][0]=RDB; // Pattern: opposite of the matching letter
2357  other_vertex_on_edge[LDB][1]=LUB;
2358  other_vertex_on_edge[LDB][2]=LDF;
2359 
2360  other_vertex_in_edge_neighbour[LDB].resize(3);
2361  other_vertex_in_edge_neighbour[LDB][0]=RUF; // Pattern: full reflection
2362  other_vertex_in_edge_neighbour[LDB][1]=RUF;
2363  other_vertex_in_edge_neighbour[LDB][2]=RUF;
2364 
2365 
2366 
2367 
2368  vertex_in_edge_neighbour[RDB].resize(3);
2369  vertex_in_edge_neighbour[RDB][0]=RUF; // Pattern: exactly one letter matches
2370  vertex_in_edge_neighbour[RDB][1]=LDF;
2371  vertex_in_edge_neighbour[RDB][2]=LUB;
2372 
2373  other_vertex_on_edge[RDB].resize(3);
2374  other_vertex_on_edge[RDB][0]=LDB; // Pattern: opposite of the matching letter
2375  other_vertex_on_edge[RDB][1]=RUB;
2376  other_vertex_on_edge[RDB][2]=RDF;
2377 
2378  other_vertex_in_edge_neighbour[RDB].resize(3);
2379  other_vertex_in_edge_neighbour[RDB][0]=LUF; // Pattern: full reflection
2380  other_vertex_in_edge_neighbour[RDB][1]=LUF;
2381  other_vertex_in_edge_neighbour[RDB][2]=LUF;
2382 
2383 
2384 
2385 
2386  vertex_in_edge_neighbour[LUB].resize(3);
2387  vertex_in_edge_neighbour[LUB][0]=LDF; // Pattern: exactly one letter matches
2388  vertex_in_edge_neighbour[LUB][1]=RUF;
2389  vertex_in_edge_neighbour[LUB][2]=RDB;
2390 
2391  other_vertex_on_edge[LUB].resize(3);
2392  other_vertex_on_edge[LUB][0]=RUB; // Pattern: opposite of the matching letter
2393  other_vertex_on_edge[LUB][1]=LDB;
2394  other_vertex_on_edge[LUB][2]=LUF;
2395 
2396  other_vertex_in_edge_neighbour[LUB].resize(3);
2397  other_vertex_in_edge_neighbour[LUB][0]=RDF; // Pattern: full reflection
2398  other_vertex_in_edge_neighbour[LUB][1]=RDF;
2399  other_vertex_in_edge_neighbour[LUB][2]=RDF;
2400 
2401 
2402 
2403 
2404  vertex_in_edge_neighbour[RUB].resize(3);
2405  vertex_in_edge_neighbour[RUB][0]=RDF; // Pattern: exactly one letter matches
2406  vertex_in_edge_neighbour[RUB][1]=LUF;
2407  vertex_in_edge_neighbour[RUB][2]=LDB;
2408 
2409  other_vertex_on_edge[RUB].resize(3);
2410  other_vertex_on_edge[RUB][0]=LUB; // Pattern: opposite of the matching letter
2411  other_vertex_on_edge[RUB][1]=RDB;
2412  other_vertex_on_edge[RUB][2]=RUF;
2413 
2414  other_vertex_in_edge_neighbour[RUB].resize(3);
2415  other_vertex_in_edge_neighbour[RUB][0]=LDF; // Pattern: full reflection
2416  other_vertex_in_edge_neighbour[RUB][1]=LDF;
2417  other_vertex_in_edge_neighbour[RUB][2]=LDF;
2418 
2419 
2420 
2421 
2422  vertex_in_edge_neighbour[LDF].resize(3);
2423  vertex_in_edge_neighbour[LDF][0]=LUB; // Pattern: exactly one letter matches
2424  vertex_in_edge_neighbour[LDF][1]=RDB;
2425  vertex_in_edge_neighbour[LDF][2]=RUF;
2426 
2427  other_vertex_on_edge[LDF].resize(3);
2428  other_vertex_on_edge[LDF][0]=RDF; // Pattern: opposite of the matching letter
2429  other_vertex_on_edge[LDF][1]=LUF;
2430  other_vertex_on_edge[LDF][2]=LDB;
2431 
2432  other_vertex_in_edge_neighbour[LDF].resize(3);
2433  other_vertex_in_edge_neighbour[LDF][0]=RUB; // Pattern: full reflection
2434  other_vertex_in_edge_neighbour[LDF][1]=RUB;
2435  other_vertex_in_edge_neighbour[LDF][2]=RUB;
2436 
2437 
2438 
2439 
2440  vertex_in_edge_neighbour[RDF].resize(3);
2441  vertex_in_edge_neighbour[RDF][0]=RUB; // Pattern: exactly one letter matches
2442  vertex_in_edge_neighbour[RDF][1]=LDB;
2443  vertex_in_edge_neighbour[RDF][2]=LUF;
2444 
2445  other_vertex_on_edge[RDF].resize(3);
2446  other_vertex_on_edge[RDF][0]=LDF; // Pattern: opposite of the matching letter
2447  other_vertex_on_edge[RDF][1]=RUF;
2448  other_vertex_on_edge[RDF][2]=RDB;
2449 
2450  other_vertex_in_edge_neighbour[RDF].resize(3);
2451  other_vertex_in_edge_neighbour[RDF][0]=LUB; // Pattern: full reflection
2452  other_vertex_in_edge_neighbour[RDF][1]=LUB;
2453  other_vertex_in_edge_neighbour[RDF][2]=LUB;
2454 
2455 
2456 
2457 
2458  vertex_in_edge_neighbour[LUF].resize(3);
2459  vertex_in_edge_neighbour[LUF][0]=LDB; // Pattern: exactly one letter matches
2460  vertex_in_edge_neighbour[LUF][1]=RUB;
2461  vertex_in_edge_neighbour[LUF][2]=RDF;
2462 
2463  other_vertex_on_edge[LUF].resize(3);
2464  other_vertex_on_edge[LUF][0]=RUF; // Pattern: opposite of the matching letter
2465  other_vertex_on_edge[LUF][1]=LDF;
2466  other_vertex_on_edge[LUF][2]=LUB;
2467 
2468  other_vertex_in_edge_neighbour[LUF].resize(3);
2469  other_vertex_in_edge_neighbour[LUF][0]=RDB; // Pattern: full reflection
2470  other_vertex_in_edge_neighbour[LUF][1]=RDB;
2471  other_vertex_in_edge_neighbour[LUF][2]=RDB;
2472 
2473 
2474 
2475 
2476  vertex_in_edge_neighbour[RUF].resize(3);
2477  vertex_in_edge_neighbour[RUF][0]=RDB; // Pattern: exactly one letter matches
2478  vertex_in_edge_neighbour[RUF][1]=LUB;
2479  vertex_in_edge_neighbour[RUF][2]=LDF;
2480 
2481  other_vertex_on_edge[RUF].resize(3);
2482  other_vertex_on_edge[RUF][0]=LUF; // Pattern: opposite of the matching letter
2483  other_vertex_on_edge[RUF][1]=RDF;
2484  other_vertex_on_edge[RUF][2]=RUB;
2485 
2486  other_vertex_in_edge_neighbour[RUF].resize(3);
2487  other_vertex_in_edge_neighbour[RUF][0]=LDB; // Pattern: full reflection
2488  other_vertex_in_edge_neighbour[RUF][1]=LDB;
2489  other_vertex_in_edge_neighbour[RUF][2]=LDB;
2490 
2491 
2492 
2493 
2494  // Loop over all vertices in the reference element
2495  for (int vertex=LDB;vertex<=RUF;vertex++)
2496  {
2497 
2498  // Loop over the three edges that are connected to this vertex
2499  for (unsigned i=0;i<3;i++)
2500  {
2501 
2502  // This is the other vertex along this edge
2503  int other_vertex=other_vertex_on_edge[vertex][i];
2504 
2505  // This is the vertex in the edge neighbour that corresponds
2506  // to the present vertex in the reference element (in
2507  // the absence of rotations)
2508  int unrotated_neigh_vertex=vertex_in_edge_neighbour[vertex][i];
2509 
2510  // This is the vertex in the edge neighbour that corresponds
2511  // to the other vertex in the reference element (in
2512  // the absence of rotations)
2513  int unrotated_neigh_other_vertex=
2514  other_vertex_in_edge_neighbour[vertex][i];
2515 
2516  // Loop over all vertices in the neighbour element
2517  for (int neigh_vertex=LDB;neigh_vertex<=RUF;neigh_vertex++)
2518  {
2519 
2520  // What rotations would turn the neigh_vertex
2521  // into the unrotated_neigh_vertex?
2522  std::set<std::pair<int,int> > vertex_rot=
2523  required_rotation[std::make_pair(neigh_vertex,
2524  unrotated_neigh_vertex)];
2525 
2526 
2527 
2528  // Loop over all "other" vertices in the neighbour element
2529  for (int neigh_other_vertex=LDB;neigh_other_vertex<=RUF;
2530  neigh_other_vertex++)
2531  {
2532 
2533  // What rotations would turn the other_neigh_vertex
2534  // into the unrotated_other_neigh_vertex?
2535  std::set<std::pair<int,int> > other_vertex_rot=
2536  required_rotation[std::make_pair(neigh_other_vertex,
2537  unrotated_neigh_other_vertex)];
2538 
2539  // What are the common rotations?
2540  std::set<std::pair<int,int> > common_rotations;
2541 
2542  // Get the intersection of the two sets
2543  std::set_intersection(vertex_rot.begin(),
2544  vertex_rot.end(),
2545  other_vertex_rot.begin(),
2546  other_vertex_rot.end(),
2547  inserter(common_rotations,
2548  common_rotations.begin()));
2549 
2550 
2551  if (common_rotations.size()>0)
2552  {
2553  for (std::set<std::pair<int,int> >::iterator it=
2554  common_rotations.begin();
2555  it!=common_rotations.end();it++)
2556  {
2557  // Copy into container
2558 
2559  // First: up equivalent:
2561  std::make_pair(std::make_pair(vertex,neigh_vertex),
2562  std::make_pair(other_vertex,
2563  neigh_other_vertex))].
2564  first=it->first;
2565 
2566  // Second: Right equivalent
2568  std::make_pair(std::make_pair(vertex,neigh_vertex),
2569  std::make_pair(other_vertex,
2570  neigh_other_vertex))].
2571  second=it->second;
2572 
2573  }
2574  }
2575 
2576  }
2577  }
2578  }
2579  }
2580 
2581 
2582 
2583 }
2584 
2585 
2586 
2587 //================================================================
2588 /// \short Is the edge neighbour (for edge "edge") specified via the pointer
2589 /// also a face neighbour for one of the two adjacent faces?
2590 //================================================================
2592  OcTree* edge_neigh_pt) const
2593 {
2594 
2595 #ifdef PARANOID
2596  // No paranoid check needed -- the default for the switch statement
2597  // catches illegal values for edge
2598 #endif
2599 
2600 
2601  // Catch stupid case: Null doesn't have a face neighbour...
2602  if (edge_neigh_pt==0)
2603  {return false;}
2604 
2605  using namespace OcTreeNames;
2606 
2607  // Auxiliary variables
2608  int face;
2609  Vector<unsigned> translate_s(3);
2610  Vector<double> s_sw(3);
2611  Vector<double> s_ne(3);
2612  int reflected_face;
2613  int diff_level;
2614 
2615  OcTree* face_neigh_pt=0;
2616 
2617  switch (edge)
2618  {
2619  case LB:
2620 
2621  // Get first face neighbour
2622  face=L;
2623  face_neigh_pt=gteq_face_neighbour(face,
2624  translate_s,
2625  s_sw,
2626  s_ne,
2627  reflected_face,
2628  diff_level);
2629 
2630  // Check if they agree...
2631  if (face_neigh_pt!=0)
2632  {
2633  if(face_neigh_pt==edge_neigh_pt)
2634  {
2635  return true;
2636  }
2637  }
2638 
2639  // Get second face neighbour
2640  face=B;
2641  face_neigh_pt=gteq_face_neighbour(face,
2642  translate_s,
2643  s_sw,
2644  s_ne,
2645  reflected_face,
2646  diff_level);
2647  // Check if they agree...
2648  if (face_neigh_pt!=0)
2649  {
2650  if (face_neigh_pt==edge_neigh_pt)
2651  {
2652  return true;
2653  }
2654  }
2655 
2656  break;
2657 
2658 
2659  case RB:
2660 
2661 
2662  // Get first face neighbour
2663  face=R;
2664  face_neigh_pt=gteq_face_neighbour(face,
2665  translate_s,
2666  s_sw,
2667  s_ne,
2668  reflected_face,
2669  diff_level);
2670 
2671  // Check if they agree...
2672  if (face_neigh_pt!=0)
2673  {
2674  if (face_neigh_pt==edge_neigh_pt)
2675  {
2676  return true;
2677  }
2678  }
2679 
2680  // Get second face neighbour
2681  face=B;
2682  face_neigh_pt=gteq_face_neighbour(face,
2683  translate_s,
2684  s_sw,
2685  s_ne,
2686  reflected_face,
2687  diff_level);
2688  // Check if they agree...
2689  if (face_neigh_pt!=0)
2690  {
2691  if (face_neigh_pt==edge_neigh_pt)
2692  {
2693  return true;
2694  }
2695  }
2696 
2697  break;
2698 
2699 
2700  case DB:
2701 
2702  // Get first face neighbour
2703  face=D;
2704  face_neigh_pt=gteq_face_neighbour(face,
2705  translate_s,
2706  s_sw,
2707  s_ne,
2708  reflected_face,
2709  diff_level);
2710 
2711  // Check if they agree...
2712  if (face_neigh_pt!=0)
2713  {
2714  if (face_neigh_pt==edge_neigh_pt)
2715  {
2716  return true;
2717  }
2718  }
2719 
2720  // Get second face neighbour
2721  face=B;
2722  face_neigh_pt=gteq_face_neighbour(face,
2723  translate_s,
2724  s_sw,
2725  s_ne,
2726  reflected_face,
2727  diff_level);
2728  // Check if they agree...
2729  if (face_neigh_pt!=0)
2730  {
2731  if (face_neigh_pt==edge_neigh_pt)
2732  {
2733  return true;
2734  }
2735  }
2736 
2737  break;
2738 
2739 
2740  case UB:
2741 
2742  // Get first face neighbour
2743  face=U;
2744  face_neigh_pt=gteq_face_neighbour(face,
2745  translate_s,
2746  s_sw,
2747  s_ne,
2748  reflected_face,
2749  diff_level);
2750 
2751  // Check if they agree...
2752  if (face_neigh_pt!=0)
2753  {
2754  if (face_neigh_pt==edge_neigh_pt)
2755  {
2756  return true;
2757  }
2758  }
2759 
2760  // Get second face neighbour
2761  face=B;
2762  face_neigh_pt=gteq_face_neighbour(face,
2763  translate_s,
2764  s_sw,
2765  s_ne,
2766  reflected_face,
2767  diff_level);
2768  // Check if they agree...
2769  if (face_neigh_pt!=0)
2770  {
2771  if (face_neigh_pt==edge_neigh_pt)
2772  {
2773  return true;
2774  }
2775  }
2776 
2777  break;
2778 
2779  case LD:
2780 
2781 
2782  // Get first face neighbour
2783  face=L;
2784  face_neigh_pt=gteq_face_neighbour(face,
2785  translate_s,
2786  s_sw,
2787  s_ne,
2788  reflected_face,
2789  diff_level);
2790 
2791  // Check if they agree...
2792  if (face_neigh_pt!=0)
2793  {
2794  if (face_neigh_pt==edge_neigh_pt)
2795  {
2796  return true;
2797  }
2798  }
2799 
2800  // Get second face neighbour
2801  face=D;
2802  face_neigh_pt=gteq_face_neighbour(face,
2803  translate_s,
2804  s_sw,
2805  s_ne,
2806  reflected_face,
2807  diff_level);
2808  // Check if they agree...
2809  if (face_neigh_pt!=0)
2810  {
2811  if (face_neigh_pt==edge_neigh_pt)
2812  {
2813  return true;
2814  }
2815  }
2816 
2817  break;
2818 
2819  case RD:
2820 
2821 
2822  // Get first face neighbour
2823  face=R;
2824  face_neigh_pt=gteq_face_neighbour(face,
2825  translate_s,
2826  s_sw,
2827  s_ne,
2828  reflected_face,
2829  diff_level);
2830 
2831  // Check if they agree...
2832  if (face_neigh_pt!=0)
2833  {
2834  if (face_neigh_pt==edge_neigh_pt)
2835  {
2836  return true;
2837  }
2838  }
2839 
2840  // Get second face neighbour
2841  face=D;
2842  face_neigh_pt=gteq_face_neighbour(face,
2843  translate_s,
2844  s_sw,
2845  s_ne,
2846  reflected_face,
2847  diff_level);
2848  // Check if they agree...
2849  if (face_neigh_pt!=0)
2850  {
2851  if (face_neigh_pt==edge_neigh_pt)
2852  {
2853  return true;
2854  }
2855  }
2856 
2857  break;
2858 
2859  case LU:
2860 
2861  // Get first face neighbour
2862  face=L;
2863  face_neigh_pt=gteq_face_neighbour(face,
2864  translate_s,
2865  s_sw,
2866  s_ne,
2867  reflected_face,
2868  diff_level);
2869 
2870  // Check if they agree...
2871  if (face_neigh_pt!=0)
2872  {
2873  if (face_neigh_pt==edge_neigh_pt)
2874  {
2875  return true;
2876  }
2877  }
2878 
2879  // Get second face neighbour
2880  face=U;
2881  face_neigh_pt=gteq_face_neighbour(face,
2882  translate_s,
2883  s_sw,
2884  s_ne,
2885  reflected_face,
2886  diff_level);
2887  // Check if they agree...
2888  if (face_neigh_pt!=0)
2889  {
2890  if (face_neigh_pt==edge_neigh_pt)
2891  {
2892  return true;
2893  }
2894  }
2895 
2896  break;
2897 
2898 
2899  case RU:
2900 
2901  // Get first face neighbour
2902  face=R;
2903  face_neigh_pt=gteq_face_neighbour(face,
2904  translate_s,
2905  s_sw,
2906  s_ne,
2907  reflected_face,
2908  diff_level);
2909 
2910  // Check if they agree...
2911  if (face_neigh_pt!=0)
2912  {
2913  if (face_neigh_pt==edge_neigh_pt)
2914  {
2915  return true;
2916  }
2917 
2918  }
2919 
2920  // Get second face neighbour
2921  face=U;
2922  face_neigh_pt=gteq_face_neighbour(face,
2923  translate_s,
2924  s_sw,
2925  s_ne,
2926  reflected_face,
2927  diff_level);
2928  // Check if they agree...
2929  if (face_neigh_pt!=0)
2930  {
2931  if (face_neigh_pt==edge_neigh_pt)
2932  {
2933  return true;
2934  }
2935  }
2936 
2937  break;
2938 
2939 
2940  case LF:
2941 
2942 
2943  // Get first face neighbour
2944  face=L;
2945  face_neigh_pt=gteq_face_neighbour(face,
2946  translate_s,
2947  s_sw,
2948  s_ne,
2949  reflected_face,
2950  diff_level);
2951 
2952  // Check if they agree...
2953  if (face_neigh_pt!=0)
2954  {
2955  if (face_neigh_pt==edge_neigh_pt)
2956  {
2957  return true;
2958  }
2959  }
2960 
2961  // Get second face neighbour
2962  face=F;
2963  face_neigh_pt=gteq_face_neighbour(face,
2964  translate_s,
2965  s_sw,
2966  s_ne,
2967  reflected_face,
2968  diff_level);
2969  // Check if they agree...
2970  if (face_neigh_pt!=0)
2971  {
2972  if (face_neigh_pt==edge_neigh_pt)
2973  {
2974  return true;
2975  }
2976  }
2977 
2978  break;
2979 
2980  case RF:
2981 
2982  // Get first face neighbour
2983  face=R;
2984  face_neigh_pt=gteq_face_neighbour(face,
2985  translate_s,
2986  s_sw,
2987  s_ne,
2988  reflected_face,
2989  diff_level);
2990 
2991  // Check if they agree...
2992  if (face_neigh_pt!=0)
2993  {
2994  if (face_neigh_pt==edge_neigh_pt)
2995  {
2996  return true;
2997  }
2998  }
2999 
3000  // Get second face neighbour
3001  face=F;
3002  face_neigh_pt=gteq_face_neighbour(face,
3003  translate_s,
3004  s_sw,
3005  s_ne,
3006  reflected_face,
3007  diff_level);
3008  // Check if they agree...
3009  if (face_neigh_pt!=0)
3010  {
3011  if (face_neigh_pt==edge_neigh_pt)
3012  {
3013  return true;
3014  }
3015  }
3016 
3017  break;
3018 
3019 
3020  case DF:
3021 
3022  // Get first face neighbour
3023  face=D;
3024  face_neigh_pt=gteq_face_neighbour(face,
3025  translate_s,
3026  s_sw,
3027  s_ne,
3028  reflected_face,
3029  diff_level);
3030 
3031  // Check if they agree...
3032  if (face_neigh_pt!=0)
3033  {
3034  if (face_neigh_pt==edge_neigh_pt)
3035  {
3036  return true;
3037  }
3038  }
3039 
3040 
3041  // Get second face neighbour
3042  face=F;
3043  face_neigh_pt=gteq_face_neighbour(face,
3044  translate_s,
3045  s_sw,
3046  s_ne,
3047  reflected_face,
3048  diff_level);
3049  // Check if they agree...
3050  if (face_neigh_pt!=0)
3051  {
3052  if (face_neigh_pt==edge_neigh_pt)
3053  {
3054  return true;
3055  }
3056  }
3057 
3058  break;
3059 
3060 
3061  case UF:
3062 
3063  // Get first face neighbour
3064  face=U;
3065  face_neigh_pt=gteq_face_neighbour(face,
3066  translate_s,
3067  s_sw,
3068  s_ne,
3069  reflected_face,
3070  diff_level);
3071 
3072  // Check if they agree...
3073  if (face_neigh_pt!=0)
3074  {
3075  if (face_neigh_pt==edge_neigh_pt)
3076  {
3077  return true;
3078  }
3079  }
3080 
3081 
3082  // Get second face neighbour
3083  face=F;
3084  face_neigh_pt=gteq_face_neighbour(face,
3085  translate_s,
3086  s_sw,
3087  s_ne,
3088  reflected_face,
3089  diff_level);
3090  // Check if they agree...
3091  if (face_neigh_pt!=0)
3092  {
3093  if (face_neigh_pt==edge_neigh_pt)
3094  {
3095  return true;
3096  }
3097  }
3098 
3099  break;
3100 
3101  default:
3102 
3103  // There is no face neighbour in this direction so they can't
3104  // agree:
3105  std::ostringstream error_stream;
3106  error_stream << "Never get here! Edge:" << Direct_string[edge]
3107  << " " << edge << std::endl;
3108  throw OomphLibError(error_stream.str(),
3109  OOMPH_CURRENT_FUNCTION,
3110  OOMPH_EXCEPTION_LOCATION);
3111 
3112 
3113  }
3114 
3115  // If we've made it to here then we've located the requested edge
3116  // but found that none of its two face neighbours match the specified
3117  // edge neighbour:
3118  return false;
3119 
3120 
3121 }
3122 
3123 
3124 //================================================================
3125 /// Find (pointer to) `greater-or-equal-sized face neighbour' in
3126 /// given direction (L/R/U/D/F/B).
3127 /// Another way of interpreting this is that we're looking for
3128 /// the neighbour across the present element's face 'direction'.
3129 /// The various arguments return additional information about the
3130 /// size and relative orientation of the neighbouring octree.
3131 /// To interpret these we use the following
3132 /// <B>General convention:</B>
3133 /// - Each face of the element that is represented by the octree
3134 /// is parametrised by two (of the three)
3135 /// local coordinates that parametrise the entire 3D element. E.g.
3136 /// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
3137 /// is parametrised by (s[0],s[2]); etc. We always identify the
3138 /// in-face coordinate with the lower (3D) index with the subscript
3139 /// _lo and the one with the larger (3D) index with the subscript _hi.
3140 /// .
3141 /// With this convention, the interpretation of the arguments is
3142 /// as follows:
3143 /// - The vector \c translate_s turns the index of the local coordinate
3144 /// in the present octree into that of the neighbour. If there are no
3145 /// rotations then \c translate_s[i] = i.
3146 /// - In the present octree, the "south west" vertex of the face
3147 /// between the present octree and its neighbour is located at
3148 /// S_lo=-1, S_hi=-1. This point is located at the (3D) local
3149 /// coordinates (\c s_sw[0], \c s_sw[1], \c s_sw[2])
3150 /// in the neighbouring octree.
3151 /// - ditto with s_ne: In the present octree, the "north east" vertex
3152 /// of the face between the present octree and its neighbour is located at
3153 /// S_lo=+1, S_hi=+1. This point is located
3154 /// at the (3D) local coordinates (\c s_ne[0], \c s_ne[1], \c s_ne[2])
3155 /// in the neighbouring octree.
3156 /// - We're looking for a neighbour in the specified \c direction. When
3157 /// viewed from the neighbouring octree, the face that separates
3158 /// the present octree from its neighbour is the neighbour's face
3159 /// \c face. If there's no rotation between the two octrees, this is a
3160 /// simple reflection: For instance, if we're looking
3161 /// for a neighhbour in the \c R [ight] \c direction, \c face will
3162 /// be \c L [eft]
3163 /// - \c diff_level <= 0 indicates the difference in refinement levels between
3164 /// the two neighbours. If \c diff_level==0, the neighbour has the
3165 /// same size as the current octree.
3166 //=================================================================
3167 OcTree* OcTree::gteq_face_neighbour(const int& direction,
3168  Vector<unsigned> &translate_s,
3169  Vector<double>& s_sw,
3170  Vector<double>& s_ne,
3171  int& face,
3172  int& diff_level) const
3173 {
3174 
3175  using namespace OcTreeNames;
3176 
3177 #ifdef PARANOID
3178  if ((direction!=L)&&(direction!=R)&&
3179  (direction!=F)&&(direction!=B)&&
3180  (direction!=U)&&(direction!=D))
3181  {
3182  std::ostringstream error_stream;
3183  error_stream << "Wrong direction: " << Direct_string[direction]
3184  << std::endl;
3185  throw OomphLibError(error_stream.str(),
3186  OOMPH_CURRENT_FUNCTION,
3187  OOMPH_EXCEPTION_LOCATION);
3188  }
3189 #endif
3190 
3191 
3192  // Maximum level to which we're allowed to descend (we only want
3193  // greater-or-equal-sized neighbours)
3194  int max_level=Level;
3195 
3196  // Current element has the following root:
3197  OcTreeRoot* orig_root_pt = dynamic_cast<OcTreeRoot*>(Root_pt);
3198 
3199  // Initialise offset in local coordinate
3200  double s_difflo=0;
3201  double s_diffhi=0;
3202 
3203  // Initialise difference in level
3204  diff_level=0;
3205 
3206  // Find neighbour
3207  OcTree* return_pt=gteq_face_neighbour(direction,s_difflo,s_diffhi,
3208  diff_level,max_level,
3209  orig_root_pt);
3210 
3211  OcTree* neighb_pt=return_pt;
3212 
3213  //Initialise the translation scheme
3214  for(unsigned i=0;i<3;i++) {translate_s[i] = i;}
3215 
3216  // If neighbour exists: What's the direction of the interfacial
3217  // face when viewed from within the neighbour element?
3218  if (neighb_pt!=0)
3219  {
3220  //Find the reflection of the original direction, which will be the
3221  //direction to the face in the neighbour, if there are no rotations.
3222  int reflected_dir = Reflect_face[direction];
3223 
3224  // These coordinates are the coordinates of the ne and sw points
3225  // in the neighbour (provided there are no rotations -- we'll correct
3226  // for these below)
3227  s_sw[0]=S_base(0,reflected_dir)+
3228  S_steplo(0,reflected_dir)*s_difflo+S_stephi(0,reflected_dir)*s_diffhi;
3229 
3230  s_sw[1]=S_base(1,reflected_dir)+
3231  S_steplo(1,reflected_dir)*s_difflo+S_stephi(1,reflected_dir)*s_diffhi;
3232 
3233  s_sw[2]=S_base(2,reflected_dir)+
3234  S_steplo(2,reflected_dir)*s_difflo+S_stephi(2,reflected_dir)*s_diffhi;
3235 
3236  s_ne[0]=S_base(0,reflected_dir)+
3237  S_steplo(0,reflected_dir)*pow(2.0,diff_level)+
3238  S_steplo(0,reflected_dir)*s_difflo+
3239  S_stephi(0,reflected_dir)*pow(2.0,diff_level)+
3240  S_stephi(0,reflected_dir)*s_diffhi;
3241 
3242  s_ne[1]=S_base(1,reflected_dir)+
3243  S_steplo(1,reflected_dir)*pow(2.0,diff_level)
3244  +S_steplo(1,reflected_dir)*s_difflo+
3245  S_stephi(1,reflected_dir)*pow(2.0,diff_level)
3246  +S_stephi(1,reflected_dir)*s_diffhi;
3247 
3248  s_ne[2]=S_base(2,reflected_dir)+
3249  S_steplo(2,reflected_dir)*pow(2.0,diff_level)
3250  +S_steplo(2,reflected_dir)*s_difflo+
3251  S_stephi(2,reflected_dir)*pow(2.0,diff_level)
3252  +S_stephi(2,reflected_dir)*s_diffhi;
3253 
3254  //If there is no rotation then the new direction is the same as the
3255  //old direction
3256  int new_dir = direction;
3257 
3258  // If necessary, rotate things around (the orientation of Up in the
3259  // neighbour might be be different from that in the present element)
3260  // If the root of the neighbour is the not same as the one of the present
3261  // element then their orientations may not be the same and the new direction
3262  // is given by :
3263  if (neighb_pt->Root_pt!=Root_pt)
3264  {
3265  new_dir=rotate(orig_root_pt->up_equivalent(
3266  neighb_pt->Root_pt),
3267  orig_root_pt->right_equivalent(
3268  neighb_pt->Root_pt),
3269  direction);
3270  }
3271 
3272  // What's the direction of the interfacial face when viewed from within
3273  // the neighbour element?
3274  face=Reflect_face[new_dir];
3275 
3276  Vector<double> s_sw_new(3), s_ne_new(3);
3277 
3278  // If the root of the present element is different from the root
3279  // of his neighbour, we have to rotate the ne and sw coordinates
3280  // to have their equivalents in the neighbour's point of view.
3281 
3282  if (neighb_pt->Root_pt!=Root_pt)
3283  {
3284  int tmp_dir;
3285  Vector<int> vect1(3);
3286  Vector<int> vect2(3);
3287  Vector<int> vect3(3);
3288  DenseMatrix<int> Mat_rot(3);
3289 
3290 
3291  // All this is just to compute the rotation matrix
3292  tmp_dir=rotate(orig_root_pt->
3293  up_equivalent(neighb_pt->Root_pt),
3294  orig_root_pt->
3295  right_equivalent(neighb_pt->Root_pt),
3296  R);
3297  vect1=Direction_to_vector[tmp_dir];
3298 
3299 
3300  tmp_dir=rotate(orig_root_pt->
3301  up_equivalent(neighb_pt->Root_pt),
3302  orig_root_pt->
3303  right_equivalent(neighb_pt->Root_pt),
3304  U);
3305  vect2=Direction_to_vector[tmp_dir];
3306 
3307 
3308  tmp_dir=rotate(orig_root_pt->
3309  up_equivalent(neighb_pt->Root_pt),
3310  orig_root_pt->
3311  right_equivalent(neighb_pt->Root_pt),
3312  F);
3313  vect3=Direction_to_vector[tmp_dir];
3314 
3315 
3316  // Setup the inverse rotation matrix
3317  for(int i=0; i<3; i++)
3318  {
3319  Mat_rot(i,0)=vect1[i];
3320  Mat_rot(i,1)=vect2[i];
3321  Mat_rot(i,2)=vect3[i];
3322  }
3323 
3324 
3325  //Initialise the translation scheme
3326  Vector<int> translate_s_new(3);
3327 
3328  // Then the rotation of the coordinates
3329  for(unsigned i=0;i<3;i++)
3330  {
3331  s_ne_new[i] = 0.0;
3332  s_sw_new[i] = 0.0;
3333  translate_s_new[i] = 0;
3334  for(unsigned k=0;k<3;k++)
3335  {
3336  s_ne_new[i] += Mat_rot(i,k)*s_ne[k];
3337  s_sw_new[i] += Mat_rot(i,k)*s_sw[k];
3338  translate_s_new[i] += Mat_rot(i,k)*translate_s[k];
3339  }
3340  }
3341 
3342  s_ne=s_ne_new;
3343  s_sw=s_sw_new;
3344 
3345  //Set the translation scheme
3346  for(unsigned i=0;i<3;i++)
3347  {
3348  // abs is ok here; not fabs!
3349  translate_s[i] = std::abs(translate_s_new[i]);
3350  }
3351  }
3352 
3353  } // end of if(neighb_pt!=0)
3354 
3355  return return_pt;
3356 
3357 }
3358 
3359 
3360 
3361 //================================================================
3362 /// Find (pointer to) `greater-or-equal-sized true edge neighbour' in
3363 /// the given direction (LB,RB,DB,UB [the back edges],
3364 /// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]).
3365 /// Another way of interpreting this is that we're looking for
3366 /// the neighbour across the present element's edge 'direction'.
3367 /// The various arguments return additional information about the
3368 /// size and relative orientation of the neighbouring octree.
3369 /// Each edge of the element that is represented by the octree
3370 /// is parametrised by one (of the three) local coordinates that
3371 /// parametrise the entire 3D element. E.g. the L[eft]B[ack] edge
3372 /// is parametrised by s[1]; the "low" vertex of this edge
3373 /// (located at the low value of this coordinate, i.e. at s[1]=-1)
3374 /// is L[eft]D[own]B[ack]. The "high" vertex of this edge (located
3375 /// at the high value of this coordinate, i.e. at s[1]=1) is
3376 /// L[eft]U[p]B[ack]; etc
3377 /// The interpretation of the arguments is as follows:
3378 /// - In a forest, an OcTree can have multiple edge neighbours
3379 /// (across an edge where multiple trees meet). \c i_root_edge_neighbour
3380 /// specifies which of these is used. Use this as "reverse communication":
3381 /// First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
3382 /// initialised to anything you want (zero, ideally). On return from
3383 /// the fct, \c n_root_edge_neighour contains the total number of true
3384 /// edge neighbours, so additional calls to the fct with
3385 /// \c i_root_edge_neighbour>0 can be made until they've all been visited.
3386 /// - The vector \c translate_s turns the index of the local coordinate
3387 /// in the present octree into that of the neighbour. If there are no
3388 /// rotations then \c translate_s[i] = i.
3389 /// - The "low" vertex of the edge in the present octree
3390 /// coincides with a certain vertex in the edge neighbour.
3391 /// In terms of the neighbour's local coordinates, this point is
3392 /// located at the (3D) local coordinates (\c s_lo[0], \c s_lo[1],
3393 /// \c s_lo[2])
3394 /// - ditto with s_hi: The "high" vertex of the edge in the present octree
3395 /// coincides with a certain vertex in the edge neighbour.
3396 /// In terms of the neighbour's local coordinates, this point is
3397 /// located at the (3D) local coordinates (\c s_hi[0], \c s_hi[1],
3398 /// \c s_hi[2])
3399 /// - We're looking for a neighbour in the specified \c direction. When
3400 /// viewed from the neighbouring octree, the edge that separates
3401 /// the present octree from its neighbour is the neighbour's edge
3402 /// \c edge. If there's no rotation between the two octrees, this is a
3403 /// simple reflection: For instance, if we're looking
3404 /// for a neighhbour in the \c DB \c direction, \c edge will
3405 /// be \c UF.
3406 /// - \c diff_level <= 0 indicates the difference in refinement levels between
3407 /// the two neighbours. If \c diff_level==0, the neighbour has the
3408 /// same size as the current octree.
3409 /// .
3410 /// \b Important: We're only looking for \b true edge neighbours
3411 /// i.e. edge neigbours that are not also face neighbours. This is an
3412 /// important difference to Samet's terminology. If the neighbour
3413 /// in a certain direction is not a true edge neighbour, or if there
3414 /// is no neighbour, then this function returns NULL.
3415 //=================================================================
3417  const unsigned& i_root_edge_neighbour,
3418  unsigned& nroot_edge_neighbour,
3419  Vector<unsigned> &translate_s,
3420  Vector<double>& s_lo,
3421  Vector<double>& s_hi,
3422  int& edge,
3423  int& diff_level) const
3424 {
3425 
3426  using namespace OcTreeNames;
3427 
3428 
3429 #ifdef PARANOID
3430  if ((direction!=LB)&&(direction!=RB)&&(direction!=DB)&&(direction!=UB)&&
3431  (direction!=LD)&&(direction!=RD)&&(direction!=LU)&&(direction!=RU)&&
3432  (direction!=LF)&&(direction!=RF)&&(direction!=DF)&&(direction!=UF))
3433  {
3434  std::ostringstream error_stream;
3435  error_stream << "Wrong direction: " << Direct_string[direction]
3436  << std::endl;
3437  throw OomphLibError(error_stream.str(),
3438  OOMPH_CURRENT_FUNCTION,
3439  OOMPH_EXCEPTION_LOCATION);
3440  }
3441 #endif
3442 
3443  // Maximum level to which we're allowed to descend (we only want
3444  // greater-or-equal-sized neighbours)
3445  int max_level=Level;
3446 
3447  // Current element has the following root:
3448  OcTreeRoot* orig_root_pt = dynamic_cast<OcTreeRoot*>(Root_pt);
3449 
3450  // Initialise offset in local coordinate along edge
3451  double s_diff=0;
3452 
3453  // Initialise difference in level
3454  diff_level=0;
3455 
3456  // Find edge neighbour
3457  OcTree* return_pt=gteq_edge_neighbour(direction,
3458  i_root_edge_neighbour,
3459  nroot_edge_neighbour,
3460  s_diff,
3461  diff_level,max_level,
3462  orig_root_pt);
3463 
3464  // Only use "true" edge neighbours
3465  if (edge_neighbour_is_face_neighbour(direction,return_pt))
3466  {
3467  return_pt=0;
3468  }
3469 
3470  // By default, we return what was returned as the true edge neighbour.
3471  OcTree* neighb_pt=return_pt;
3472 
3473  //Initialise the translation scheme
3474  for(unsigned i=0;i<3;i++) {translate_s[i] = i;}
3475 
3476  // If neighbour exists: What's the direction of the interfacial
3477  // edge when viewed from within the neighbour element?
3478  if (neighb_pt!=0)
3479  {
3480  //Find the reflection of the original direction, which will be the
3481  //direction to the edge in the neighbour, if there are no rotations.
3482  int reflected_dir = Reflect_edge[direction];
3483 
3484  // These coordinates are the coordinates of the "low" and "high" points
3485  // in the neighbour (provided there are no rotations -- we'll correct
3486  // for these below)
3487  s_lo[0]=S_base_edge(0,reflected_dir)+
3488  S_step_edge(0,reflected_dir)*s_diff;
3489 
3490  s_lo[1]=S_base_edge(1,reflected_dir)+
3491  S_step_edge(1,reflected_dir)*s_diff;
3492 
3493  s_lo[2]=S_base_edge(2,reflected_dir)+
3494  S_step_edge(2,reflected_dir)*s_diff;
3495 
3496  s_hi[0]=S_base_edge(0,reflected_dir)+
3497  S_step_edge(0,reflected_dir)*pow(2.0,diff_level)+
3498  S_step_edge(0,reflected_dir)*s_diff;
3499 
3500  s_hi[1]=S_base_edge(1,reflected_dir)+
3501  S_step_edge(1,reflected_dir)*pow(2.0,diff_level)
3502  +S_step_edge(1,reflected_dir)*s_diff;
3503 
3504  s_hi[2]=S_base_edge(2,reflected_dir)+
3505  S_step_edge(2,reflected_dir)*pow(2.0,diff_level)
3506  +S_step_edge(2,reflected_dir)*s_diff;
3507 
3508  //If there is no rotation then the new direction is the same as the
3509  //old direction
3510  int new_dir = direction;
3511 
3512 
3513  // If necessary, rotate things around (the orientation of Up in the
3514  // neighbour might be be different from that in the present element)
3515  // If the root of the neighbour is the not same as the one of the present
3516  // element then their orientations may not be the same and the new
3517  // direction is given by :
3518  if (neighb_pt->Root_pt!=Root_pt)
3519  {
3520  int new_up=orig_root_pt->
3521  up_equivalent(neighb_pt->Root_pt);
3522 
3523  int new_right=orig_root_pt->
3524  right_equivalent(neighb_pt->Root_pt);
3525 
3526  new_dir=rotate(new_up,new_right,direction);
3527  }
3528 
3529  // What's the direction of the interfacial edge when viewed from within
3530  // the neighbour element (including rotations!)
3531  edge=Reflect_edge[new_dir];
3532 
3533  // Get ready to rotate the local coordinates
3534  Vector<double> s_lo_new(3), s_hi_new(3);
3535 
3536  // If the root of the present element is different from the root
3537  // of his neighbour, we have to rotate the lo and hi coordinates
3538  // to have their equivalents from the neighbour's point of view.
3539  if ((neighb_pt->Root_pt!=Root_pt))
3540  {
3541 
3542  int tmp_dir;
3543  Vector<int> vect1(3);
3544  Vector<int> vect2(3);
3545  Vector<int> vect3(3);
3546  DenseMatrix<int> Mat_rot(3);
3547 
3548  // All this is just to compute the rotation matrix
3549  tmp_dir = rotate(orig_root_pt->up_equivalent(
3550  neighb_pt->Root_pt),
3551  orig_root_pt->right_equivalent(
3552  neighb_pt->Root_pt),
3553  R);
3554  vect1=Direction_to_vector[tmp_dir];
3555 
3556 
3557  tmp_dir = rotate(orig_root_pt->up_equivalent(
3558  neighb_pt->Root_pt),
3559  orig_root_pt->right_equivalent(
3560  neighb_pt->Root_pt),
3561  U);
3562  vect2=Direction_to_vector[tmp_dir];
3563 
3564 
3565  tmp_dir = rotate(orig_root_pt->up_equivalent(
3566  neighb_pt->Root_pt),
3567  orig_root_pt->right_equivalent(
3568  neighb_pt->Root_pt),
3569  F);
3570  vect3=Direction_to_vector[tmp_dir];
3571 
3572 
3573  // Setup the inverse rotation matrix
3574  for(int i=0; i<3; i++)
3575  {
3576  Mat_rot(i,0)=vect1[i];
3577  Mat_rot(i,1)=vect2[i];
3578  Mat_rot(i,2)=vect3[i];
3579  }
3580 
3581  //Initialise the translation scheme
3582  Vector<int> translate_s_new(3);
3583 
3584  // Then the rotation of the coordinates
3585  for(unsigned i=0;i<3;i++)
3586  {
3587  s_hi_new[i] = 0.0;
3588  s_lo_new[i] = 0.0;
3589  translate_s_new[i] = 0;
3590  for(unsigned k=0;k<3;k++)
3591  {
3592  s_hi_new[i] += Mat_rot(i,k)*s_hi[k];
3593  s_lo_new[i] += Mat_rot(i,k)*s_lo[k];
3594  translate_s_new[i] += Mat_rot(i,k)*translate_s[k];
3595  }
3596  }
3597 
3598  s_lo=s_lo_new;
3599  s_hi=s_hi_new;
3600 
3601  //Set the translation scheme
3602  for(unsigned i=0;i<3;i++)
3603  {
3604  // abs is ok here; not fabs!
3605  translate_s[i] = std::abs(translate_s_new[i]);
3606  }
3607  }
3608 
3609  } // end if for (neighb_pt!=0)
3610 
3611  return return_pt;
3612 
3613 }
3614 
3615 
3616 
3617 
3618 //==========================================================================
3619 /// Find `greater-or-equal-sized face neighbour' in given direction
3620 /// (L/R/U/D/B/F).
3621 ///
3622 /// This is an auxiliary routine which allows neighbour finding in adjacent
3623 /// octrees. Needs to keep track of previous son types and
3624 /// the maximum level to which search is performed.
3625 ///
3626 /// Parameters:
3627 ///
3628 /// - direction: (L/R/U/D/B/F) Direction in which neighbour has to be found.
3629 /// - s_difflo/s_diffhi: Offset of left/down/back vertex from
3630 /// corresponding vertex in
3631 /// neighbour. Note that this is input/output as it needs to be incremented/
3632 /// decremented during the recursive calls to this function.
3633 /// - face: We're looking for the neighbour across our face 'direction'
3634 /// (L/R/U/D/B/F). When viewed from the neighbour, this face is
3635 /// `face' (L/R/U/D/B/F). [If there's no relative rotation between neighbours
3636 /// then this is a mere reflection, e.g. direction=F --> face=B etc.]
3637 /// - diff_level <= 0 indicates the difference in octree levels
3638 /// between the current element and its neighbour.
3639 /// - max_level is the maximum level to which the neighbour search is
3640 /// allowed to proceed. This is necessary because in a forest,
3641 /// the neighbour search isn't based on pure recursion.
3642 /// - orig_root_pt identifies the root node of the element whose
3643 /// neighbour we're really trying to find by all these recursive calls.
3644 //===========================================================================
3645 OcTree* OcTree::gteq_face_neighbour(const int& direction,
3646  double& s_difflo,
3647  double& s_diffhi,
3648  int& diff_level,
3649  int max_level,
3650  OcTreeRoot* orig_root_pt) const
3651 {
3652 
3653 using namespace OcTreeNames;
3654 
3655 
3656 #ifdef PARANOID
3657  if ((direction!=L)&&(direction!=R)&&
3658  (direction!=F)&&(direction!=B)&&
3659  (direction!=U)&&(direction!=D))
3660  {
3661  std::ostringstream error_stream;
3662  error_stream << "Wrong direction: " << Direct_string[direction]
3663  << std::endl;
3664  throw OomphLibError(error_stream.str(),
3665  OOMPH_CURRENT_FUNCTION,
3666  OOMPH_EXCEPTION_LOCATION);
3667  }
3668 #endif
3669 
3670  OcTree* next_el_pt;
3671  OcTree* return_el_pt;
3672 
3673 
3674  // STEP 1: Find the neighbour's father
3675  //--------
3676 
3677  // Does the element have a father?
3678  if (Father_pt!=0)
3679  {
3680 
3681  // If the present octant (whose location inside its
3682  // father element is specified by Son_type) is adjacent to the
3683  // father's face in the required direction, then its neighbour has
3684  // a different father ---> we need to climb up the tree to
3685  // the father and find his neighbour in the required direction
3686  if (Is_adjacent(direction,Son_type))
3687  {
3688 
3689  next_el_pt=dynamic_cast<OcTree*>(Father_pt)->
3690  gteq_face_neighbour(direction,
3691  s_difflo, s_diffhi,diff_level,max_level,
3692  orig_root_pt);
3693  }
3694 
3695 
3696  // If the present octant is not adjacent to the
3697  // father's face in the required direction, then the
3698  // neighbour has the same father and will be obtained
3699  // by the appropriate reflection inside the father element
3700  else
3701  {
3702  next_el_pt = dynamic_cast<OcTree*>(Father_pt);
3703  }
3704 
3705 
3706  // We're about to ascend one level:
3707  diff_level-=1;
3708 
3709  // Work out position of lower-left corner of present face
3710  // in its father element
3711 
3712  s_difflo+=pow(0.5,-diff_level)*S_directlo(direction,Son_type);
3713  s_diffhi+=pow(0.5,-diff_level)*S_directhi(direction,Son_type);
3714 
3715  // STEP 2: We have now located the neighbour's father and need to
3716  // ------- find the appropriate son.
3717 
3718  // Buffer cases where the neighbour (and hence its father) lie outside
3719  // the boundary
3720  if (next_el_pt!=0)
3721  {
3722 
3723  // If the father is a leaf then we can't descend to the same
3724  // level as the present node ---> simply return the father himself
3725  // as the (greater) neighbour. Same applies if we are about
3726  // to descend lower than the max_level (in a neighbouring tree)
3727  if ((next_el_pt->Son_pt.size()==0)||(next_el_pt->Level>max_level-1))
3728  {
3729  return_el_pt=next_el_pt;
3730  }
3731 
3732  // We have located the neighbour's father: The position of the
3733  // neighbour is obtained by `reflecting' the position of the
3734  // node itself.
3735  else
3736  {
3737 
3738  // By default (in the absence of rotations) we obtain the
3739  // son octant by reflecting
3740  int son_octant=Reflect(direction,Son_type);
3741 
3742 
3743  // If there is a rotation, we rotate the son octant
3744  if (orig_root_pt!=next_el_pt->Root_pt)
3745  {
3746  int my_up=dynamic_cast<OcTreeRoot*>(Root_pt)
3747  ->up_equivalent(next_el_pt->Root_pt);
3748 
3749  int my_right=dynamic_cast<OcTreeRoot*>(Root_pt)
3750  ->right_equivalent(next_el_pt->Root_pt);
3751 
3752  son_octant=rotate(my_up,my_right,son_octant);
3753  }
3754 
3755  return_el_pt=dynamic_cast<OcTree*>(next_el_pt->Son_pt[son_octant]);
3756 
3757 
3758  // Work out position of lower-left corner of present face
3759  // in next higher element
3760  s_difflo-=pow(0.5,-diff_level)*S_directlo(direction,Son_type);
3761  s_diffhi-=pow(0.5,-diff_level)*S_directhi(direction,Son_type);
3762 
3763  // We have just descended one level
3764  diff_level+=1;
3765  }
3766  }
3767  // The neighbour's father lies outside the boundary --> the neighbour
3768  // itself does too --> return NULL.
3769  else
3770  {
3771  return_el_pt=0;
3772  }
3773 
3774  }
3775 
3776  // Element does not have a father --> check if it has a neighbouring
3777  // tree in the appropriate direction
3778  else
3779  {
3780  // Find neighbouring root
3781  if(Root_pt->neighbour_pt(direction)!=0)
3782  {
3783  return_el_pt=dynamic_cast<OcTree*>(Root_pt->neighbour_pt(direction));
3784  }
3785 
3786  //No neighbouring tree, so there really is no neighbour --> return NULL
3787  else
3788  {
3789  return_el_pt=0;
3790  }
3791  }
3792 
3793  return return_el_pt;
3794 
3795 }
3796 
3797 
3798 
3799 
3800 //==========================================================================
3801 /// Find `greater-or-equal-sized edge neighbour' in given direction
3802 /// (LB,RB,DB,UB [the back edges],
3803 /// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]).
3804 ///
3805 /// This is an auxiliary routine which allows neighbour finding in adjacent
3806 /// octrees. Needs to keep track of previous son types and
3807 /// the maximum level to which search is performed.
3808 ///
3809 /// Parameters:
3810 ///
3811 /// - direction: (LB,RB/...) Direction in which neighbour has to be found.
3812 /// - In a forest, an OcTree can have multiple edge neighbours
3813 /// (across an edge where multiple trees meet). \c i_root_edge_neighbour
3814 /// specifies which of these is used. Use this as "reverse communication":
3815 /// First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
3816 /// initialised to anything you want (zero, ideally). On return from
3817 /// the fct, \c n_root_edge_neighour contains the total number of true
3818 /// edge neighbours, so additional calls to the fct with
3819 /// \c i_root_edge_neighbour>0 can be made until they've all been visited.
3820 /// - s_diff: Offset of left/down/back vertex from
3821 /// corresponding vertex in
3822 /// neighbour. Note that this is input/output as it needs to be incremented/
3823 /// decremented during the recursive calls to this function.
3824 /// - diff_level <= 0 indicates the difference in octree levels
3825 /// between the current element and its neighbour.
3826 /// - max_level is the maximum level to which the neighbour search is
3827 /// allowed to proceed. This is necessary because in a forest,
3828 /// the neighbour search isn't based on pure recursion.
3829 /// - orig_root_pt identifies the root node of the element whose
3830 /// neighbour we're really trying to find by all these recursive calls.
3831 //===========================================================================
3832 OcTree* OcTree::gteq_edge_neighbour(const int& direction,
3833  const unsigned& i_root_edge_neighbour,
3834  unsigned& nroot_edge_neighbour,
3835  double& s_diff,
3836  int& diff_level,
3837  int max_level,
3838  OcTreeRoot* orig_root_pt) const
3839 {
3840 
3841 using namespace OcTreeNames;
3842 
3843 
3844 #ifdef PARANOID
3845  if ((direction!=LB)&&(direction!=RB)&&(direction!=DB)&&(direction!=UB)&&
3846  (direction!=LD)&&(direction!=RD)&&(direction!=LU)&&(direction!=RU)&&
3847  (direction!=LF)&&(direction!=RF)&&(direction!=DF)&&(direction!=UF))
3848  {
3849  std::ostringstream error_stream;
3850  error_stream << "Wrong direction: " << Direct_string[direction]
3851  << std::endl;
3852  throw OomphLibError(error_stream.str(),
3853  OOMPH_CURRENT_FUNCTION,
3854  OOMPH_EXCEPTION_LOCATION);
3855  }
3856 #endif
3857 
3858  // Initialise total number of edge neighbours available across
3859  // edges of octree roots
3860  nroot_edge_neighbour=0;
3861 
3862 
3863  OcTree* next_el_pt=0;
3864  OcTree* return_el_pt=0;
3865 
3866 
3867  // STEP 1: Find the common ancestor
3868  //--------
3869 
3870  // Does the element have a father?
3871  if (Father_pt!=0)
3872  {
3873 
3874 
3875  // If the present octant (whose location inside its
3876  // father element is specified by Son_type) is adjacent to the
3877  // father's edge in the required direction, then its neighbour has
3878  // a different father ---> we need to climb up the tree to
3879  // the father and find his neighbour in the required direction
3880  if (Is_adjacent(direction,Son_type))
3881  {
3882  next_el_pt=dynamic_cast<OcTree*>(Father_pt)->
3883  gteq_edge_neighbour(direction,
3884  i_root_edge_neighbour,
3885  nroot_edge_neighbour,
3886  s_diff,diff_level,max_level,
3887  orig_root_pt);
3888  }
3889  // If the present octant (whose location inside its
3890  // father element is specified by Son_type) is adjacent to the
3891  // father's face in the required direction, then its neighbour has
3892  // a different father ---> we need to climb up the tree to
3893  // the father and find his neighbour in the required direction,
3894  // crossing the face as we do so.
3895  else if (Common_face(direction,Son_type)!=OMEGA)
3896  {
3897  // We're going across a face:
3898  double s_difflo=0.0;
3899  double s_diffhi=0.0;
3900  int diff_level_edge=0;
3901 
3902  next_el_pt=dynamic_cast<OcTree*>(Father_pt)->
3904  s_difflo, s_diffhi,
3905  diff_level_edge,max_level,
3906  orig_root_pt);
3907  }
3908 
3909  // If the present octant is not adjacent to the
3910  // father's face/edge in the required direction, then the
3911  // neighbour has the same father and will be obtained
3912  // by the appropriate reflection inside the father element
3913  else
3914  {
3915  next_el_pt = dynamic_cast<OcTree*>(Father_pt);
3916  }
3917 
3918 
3919  // We're about to ascend one level:
3920  diff_level-=1;
3921 
3922  // Work out position of "low" vertex of present edge
3923  // in its father element
3924  s_diff+=pow(0.5,-diff_level)*S_direct_edge(direction,Son_type);
3925 
3926  // STEP 2: We have now located the neighbour's father and need to
3927  // ------- find the appropriate son.
3928 
3929  // Buffer cases where ...
3930  if (next_el_pt!=0)
3931  {
3932 
3933  // If the father is a leaf then we can't descend to the same
3934  // level as the present node ---> simply return the father himself
3935  // as the (greater) neighbour. Same applies if we are about
3936  // to descend lower than the max_level (in a neighbouring tree)
3937  if ((next_el_pt->Son_pt.size()==0)||(next_el_pt->Level>max_level-1))
3938  {
3939  return_el_pt=next_el_pt;
3940  }
3941 
3942  // We have located the neighbour's father: The position of the
3943  // neighbour is obtained by `reflecting' the position of the
3944  // node itself.
3945  else
3946  {
3947 
3948  // By default (in the absence of rotations) we obtain the
3949  // son octant by reflecting
3950  int son_octant=Reflect(direction,Son_type);
3951 
3952  // If there is a rotation, we rotate the son octant
3953  if (orig_root_pt!=next_el_pt->Root_pt)
3954  {
3955  int my_up=dynamic_cast<OcTreeRoot*>(Root_pt)
3956  ->up_equivalent(next_el_pt->Root_pt);
3957 
3958  int my_right=dynamic_cast<OcTreeRoot*>(Root_pt)
3959  ->right_equivalent(next_el_pt->Root_pt);
3960 
3961  son_octant=rotate(my_up,my_right,son_octant);
3962  }
3963 
3964  return_el_pt=dynamic_cast<OcTree*>(next_el_pt->Son_pt[son_octant]);
3965 
3966  // Work out position of "low" vertex of present edge
3967  // in next higher element
3968  s_diff-=pow(0.5,-diff_level)*S_direct_edge(direction,Son_type);
3969 
3970  // We have just descended one level
3971  diff_level+=1;
3972  }
3973  }
3974  // The neighbour's father lies outside the boundary --> the neighbour
3975  // itself does too --> return NULL.
3976  else
3977  {
3978  return_el_pt=0;
3979  }
3980 
3981  }
3982 
3983  // Element does not have a father --> check if it has a neighbouring
3984  // tree in the appropriate direction
3985  else
3986  {
3987  // Get total number of edge neighbours available across
3988  // edges of octree roots for return
3989  nroot_edge_neighbour=
3990  dynamic_cast<OcTreeRoot*>(Root_pt)->nedge_neighbour(direction);
3991 
3992  // Get vector of edge neighbours (if any) in appropriate direction
3993  Vector<TreeRoot*> edge_neighbour_pt=
3994  dynamic_cast<OcTreeRoot*>(Root_pt)->edge_neighbour_pt(direction);
3995  unsigned n_neigh=edge_neighbour_pt.size();
3996  if (n_neigh>0)
3997  {
3998  return_el_pt=dynamic_cast<OcTree*>(
3999  edge_neighbour_pt[i_root_edge_neighbour]);
4000  }
4001  else
4002  {
4003  return_el_pt=0;
4004  }
4005  }
4006 
4007  return return_el_pt;
4008 
4009 }
4010 
4011 
4012 
4013 
4014 //================================================================
4015 /// Self-test: Check neighbour finding routine. For each element
4016 /// in the tree and for each vertex, determine the
4017 /// distance between the vertex and its position in the
4018 /// neigbour. . If the difference is less than
4019 /// Tree::Max_neighbour_finding_tolerance.
4020 /// return success (0), otherwise failure (1)
4021 //=================================================================
4023 {
4024 
4025  // Stick pointers to all nodes into Vector and number elements
4026  // in the process
4027  Vector<Tree*> all_nodes_pt;
4028  stick_all_tree_nodes_into_vector(all_nodes_pt);
4029 
4030  long int count=0;
4031  unsigned long num_nodes=all_nodes_pt.size();
4032 
4033  for (unsigned long i=0;i<num_nodes;i++)
4034  {
4035  all_nodes_pt[i]->object_pt()->set_number(++count);
4036  }
4037 
4038  // Check neighbours vertices and their opposite points
4039  std::ofstream neighbours_file;
4040  std::ofstream no_true_edge_file;
4041  std::ofstream neighbours_txt_file;
4042 
4043  double max_error_face=0.0;
4044  OcTree::doc_face_neighbours(all_nodes_pt,neighbours_file,
4045  neighbours_txt_file, max_error_face);
4046 
4047  double max_error_edge=0.0;
4048  OcTree::doc_true_edge_neighbours(all_nodes_pt,neighbours_file,
4049  no_true_edge_file,
4050  neighbours_txt_file, max_error_edge);
4051 
4052 
4053  bool failed=false;
4054  if (max_error_face>Max_neighbour_finding_tolerance)
4055  {
4056  oomph_info
4057  << "\n \n Failed self_test() for OcTree because of faces: Max. error "
4058  << max_error_face << std::endl<< std::endl;
4059  failed=true;
4060  }
4061 
4062  if (max_error_edge>Max_neighbour_finding_tolerance)
4063  {
4064  oomph_info
4065  << "\n \n Failed self_test() for OcTree because of edges: Max. error "
4066  << max_error_edge << std::endl<< std::endl;
4067  failed=true;
4068  }
4069 
4070  double max_error=max_error_face;
4071  if (max_error_edge>max_error) max_error=max_error_edge;
4072 
4073  if (failed)
4074  {
4075  return 1;
4076  }
4077  else
4078  {
4079  oomph_info << "Passed self_test() for OcTree: Max. error "
4080  << max_error << std::endl;
4081  return 0;
4082  }
4083 
4084 }
4085 
4086 
4087 
4088 
4089 //=================================================================
4090 /// Doc/check all face neighbours of octree (nodes) contained in the
4091 /// Vector forest_node_pt. Output into neighbours_file which can
4092 /// be viewed from tecplot with OcTreeNeighbours.mcr
4093 /// Neighbour info and errors are displayed on
4094 /// neighbours_txt_file. Finally, compute the max. error between
4095 /// vertices when viewed from neighhbouring element.
4096 /// If the two filestreams are closed, output is suppressed.
4097 /// (Static function.)
4098 //=================================================================
4100  std::ofstream& neighbours_file,
4101  std::ofstream& neighbours_txt_file,
4102  double& max_error)
4103 {
4104  using namespace OcTreeNames;
4105 
4106  int diff_level;
4107  int face=OMEGA;
4108 
4109  Vector<double> s(3);
4110  Vector<double> x(3);
4111 
4112  Vector<double> s_sw(3);
4113  Vector<double> s_ne(3);
4114  Vector<unsigned> translate_s(3);
4115 
4116  Vector<double> x_small(3);
4117  Vector<double> x_large(3);
4118 
4119 
4120  // Initialise error in vertex positions
4121  max_error=0.0;
4122 
4123  // Loop over all elements
4124  // ----------------------
4125  unsigned long num_nodes=forest_nodes_pt.size();
4126 
4127  for (unsigned long i=0;i<num_nodes;i++)
4128  {
4129 
4130 
4131  // Doc the element itself
4132  OcTree* el_pt = dynamic_cast<OcTree*>(forest_nodes_pt[i]);
4133 
4134  //If the object is incomplete omit
4135  if(el_pt->object_pt()->nodes_built())
4136  {
4137  // Print it
4138  if (neighbours_file.is_open())
4139  {
4140  neighbours_file << "#---------------------------------" << std::endl;
4141  neighbours_file << "#The element itself: " << i << std::endl;
4142  neighbours_file << "#---------------------------------" << std::endl;
4143  dynamic_cast<RefineableQElement<3>*>(
4144  el_pt->object_pt())->output_corners(neighbours_file,"BLACK");
4145  }
4146 
4147  // Loop over directions to find neighbours
4148  // ----------------------------------------
4149  for (int direction=L;direction<=F;direction++)
4150  {
4151 
4152  // Initialise difference in levels and coordinate offset
4153  diff_level=0;
4154 
4155  // Find greater-or-equal-sized neighbour...
4156  OcTree* neighb_pt=el_pt->gteq_face_neighbour(direction, translate_s,
4157  s_sw, s_ne,
4158  face,diff_level);
4159 
4160  // If neighbour exists and nodes are created
4161  if((neighb_pt!=0) && (neighb_pt->object_pt()->nodes_built()))
4162  {
4163 
4164  // Doc neighbour stats
4165  if (neighbours_txt_file.is_open())
4166  {
4167  neighbours_txt_file<< Direct_string[direction]
4168  << " neighbour of "
4169  << el_pt->object_pt()->number() << " is "
4170  << neighb_pt->object_pt()->number()
4171  << " diff_level " << diff_level
4172  << ". Inside the neighbour the face is "
4173  << Direct_string[face] << std::endl;
4174  }
4175 
4176  // Plot neighbour in the appropriate colour
4177  if (neighbours_file.is_open())
4178  {
4179  neighbours_file << "#---------------------------------" << std::endl;
4180  neighbours_file << "#Neighbour element: " << Direct_string[direction]
4181  << std::endl;
4182  neighbours_file << "#---------------------------------" << std::endl;
4183  dynamic_cast<RefineableQElement<3>*>(
4184  neighb_pt->object_pt())->
4185  output_corners(neighbours_file,Colour[direction]);
4186  }
4187 
4188  // Check that local coordinates in the larger element
4189  // lead to the same spatial point as the node vertices
4190  // in the current element
4191  if (neighbours_file.is_open())
4192  {
4193  neighbours_file << "ZONE I=2 C=" << Colour[direction] << std::endl;
4194  }
4195 
4196  // "South west" vertex in the interfacial face
4197  //--------------------------------------------
4198 
4199  // Get coordinates in large (neighbour) element
4200  s[0]=s_sw[0];
4201  s[1]=s_sw[1];
4202  s[2]=s_sw[2];
4203  neighb_pt->object_pt()->get_x(s,x_large);
4204 
4205  // Get coordinates in small element
4206  Vector<double> s(3);
4207  s[0]=S_base(0,direction);
4208  s[1]=S_base(1,direction);
4209  s[2]=S_base(2,direction);
4210  el_pt->object_pt()->get_x(s,x_small);
4211 
4212  double error=
4213  sqrt(pow(x_small[0]-x_large[0],2)+
4214  pow(x_small[1]-x_large[1],2)+
4215  pow(x_small[2]-x_large[2],2));
4216 
4217  if (neighbours_txt_file.is_open())
4218  {
4219  neighbours_txt_file << "Error (1) " << error << std::endl;
4220  }
4221 
4222  if (std::fabs(error)>max_error)
4223  {
4224  max_error=std::fabs(error);
4225  }
4226 
4227 
4228  // Check error and doc mismatch if required
4229  bool stop=false;
4230  std::ofstream mismatch_file;
4231  if (std::fabs(error)>Max_neighbour_finding_tolerance)
4232  {
4233  stop=true;
4234  mismatch_file.open("mismatch.dat");
4235  mismatch_file << "ZONE" << std::endl;
4236  mismatch_file<< x_large[0]<<" "<<x_large[1]<<" "<<x_large[2]<<" 2 \n";
4237  mismatch_file<< x_small[0]<<" "<<x_small[1]<<" "<<x_small[2]<<" 3 \n";
4238  }
4239 
4240  if (neighbours_file.is_open())
4241  {
4242  neighbours_file << "#SOUTH WEST: " << std::endl;
4243  neighbours_file << x_large[0] << " "
4244  << x_large[1] << " "
4245  << x_large[2] <<" 40 \n";
4246  }
4247 
4248 
4249  // "North east" vertex in the interfacial face
4250  //--------------------------------------------
4251 
4252  // Get coordinates in large (neighbour) element
4253  s[0]=s_ne[0];
4254  s[1]=s_ne[1];
4255  s[2]=s_ne[2];
4256  neighb_pt->object_pt()->get_x(s,x_large);
4257 
4258  //Get coordinates in small element
4259  s[0]=S_base(0,direction)+S_steplo(0,direction)+S_stephi(0,direction);
4260  s[1]=S_base(1,direction)+S_steplo(1,direction)+S_stephi(1,direction);
4261  s[2]=S_base(2,direction)+S_steplo(2,direction)+S_stephi(2,direction);
4262  el_pt->object_pt()->get_x(s,x_small);
4263 
4264  // Output
4265  if (neighbours_file.is_open())
4266  {
4267  neighbours_file << "#NORTH EAST: " << std::endl;
4268  neighbours_file << x_large[0]
4269  << " " << x_large[1]
4270  << " " << x_large[2]
4271  << " 80 \n";
4272  }
4273 
4274 
4275  error=
4276  sqrt(pow(x_small[0]-x_large[0],2)+
4277  pow(x_small[1]-x_large[1],2)+
4278  pow(x_small[2]-x_large[2],2));
4279  if (neighbours_txt_file.is_open())
4280  {
4281  neighbours_txt_file << "Error (2) " << error << std::endl;
4282  }
4283 
4284  if (std::fabs(error)>max_error)
4285  {
4286  max_error=std::fabs(error);
4287  }
4288 
4289 
4290  // Check error and doc mismatch if required
4291  if (std::fabs(error)>Max_neighbour_finding_tolerance) stop=true;
4292  if (stop)
4293  {
4294  if (!mismatch_file.is_open())
4295  {
4296  mismatch_file.open("mismatch.dat");
4297  }
4298  mismatch_file << "ZONE" << std::endl;
4299  mismatch_file << x_large[0] << " " << x_large[1] << " "
4300  << x_large[2] <<" 2 \n";
4301  mismatch_file << x_small[0] << " " << x_small[1] << " "
4302  << x_small[2] <<" 3 \n";
4303  mismatch_file.close();
4304  oomph_info << "Error for direction " << Direct_string[direction]
4305  << " with element "<<i+1<< std::endl;
4306  pause("Error");
4307  }
4308  }
4309 
4310  // If neighbour does not exist: Insert blank zones into file
4311  // so that tecplot can find six neighbours for every element
4312  else
4313  {
4314  if (neighbours_file.is_open())
4315  {
4316 
4317  neighbours_file << "#---------------------------------" << std::endl;
4318  neighbours_file << "#no Neighbour in direct: "
4319  << Direct_string[direction]
4320  << std::endl;
4321  neighbours_file << "#---------------------------------" << std::endl;
4322 
4323  dynamic_cast<RefineableQElement<3>*>(
4324  el_pt->object_pt())->output_corners(neighbours_file,"WHITE");
4325  neighbours_file << "ZONE I=2 \n";
4326  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4327  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4328  }
4329  }
4330 
4331  if (neighbours_file.is_open())
4332  {
4333  neighbours_file << std::endl << std::endl << std::endl;
4334  }
4335  }
4336  } //End of case when element can be documented
4337  }
4338 }
4339 
4340 
4341 //////////////////////////////////////////////////////////////////////
4342 //////////////////////////////////////////////////////////////////////
4343 
4344 
4345 
4346 //=================================================================
4347 /// Doc/check all true edge neighbours of octree (nodes) contained
4348 /// in the Vector forest_node_pt. Output into neighbours_file which can
4349 /// be viewed from tecplot with OcTreeNeighbours.mcr
4350 /// Neighbour info and errors are displayed on
4351 /// neighbours_txt_file. Finally, compute the max. error between
4352 /// vertices when viewed from neighhbouring element.
4353 /// If the two filestreams are closed, output is suppressed.
4354 /// (Static function).
4355 //=================================================================
4357  std::ofstream& neighbours_file,
4358  std::ofstream& no_true_edge_file,
4359  std::ofstream& neighbours_txt_file,
4360  double& max_error)
4361 {
4362  using namespace OcTreeNames;
4363 
4364  int diff_level;
4365  int edge=OMEGA;
4366 
4367  Vector<double> s(3);
4368  Vector<double> x(3);
4369 
4370  Vector<double> s_lo(3);
4371  Vector<double> s_hi(3);
4372  Vector<unsigned> translate_s(3);
4373 
4374  Vector<double> x_small(3);
4375  Vector<double> x_large(3);
4376 
4377 
4378  // Initialise error in vertex positions
4379  max_error=0.0;
4380 
4381  // Loop over all elements
4382  // ----------------------
4383  unsigned long num_nodes=forest_nodes_pt.size();
4384 
4385  for (unsigned long i=0;i<num_nodes;i++)
4386  {
4387 
4388  // Doc the element itself
4389  OcTree* el_pt = dynamic_cast<OcTree*>(forest_nodes_pt[i]);
4390 
4391  //If the object is incomplete omit
4392  if(el_pt->object_pt()->nodes_built())
4393  {
4394  // Print it
4395  if (neighbours_file.is_open())
4396  {
4397  neighbours_file << "#---------------------------------" << std::endl;
4398  neighbours_file << "#The element itself: " << i << std::endl;
4399  neighbours_file << "#---------------------------------" << std::endl;
4400  dynamic_cast<RefineableQElement<3>*>(
4401  el_pt->object_pt())->output_corners(neighbours_file,"BLACK");
4402  }
4403 
4404  // Loop over directions to find edge neighbours
4405  // --------------------------------------------
4406  for (int direction=LB;direction<=UF;direction++)
4407  {
4408  // Initialise difference in levels
4409  diff_level=0;
4410 
4411  // For now simply doc the zero-th edge neighbour (if any)
4412  unsigned i_root_edge_neighbour=0;
4413  unsigned nroot_edge_neighbour=0;
4414 
4415  // Find greater-or-equal-sized edge neighbour...
4416  OcTree* neighb_pt=el_pt->gteq_true_edge_neighbour(direction,
4417  i_root_edge_neighbour,
4418  nroot_edge_neighbour,
4419  translate_s,
4420  s_lo, s_hi,
4421  edge,diff_level);
4422 
4423  // If neighbour exists and nodes are created
4424  if((neighb_pt!=0) && (neighb_pt->object_pt()->nodes_built()))
4425  {
4426 
4427  // Doc neighbour stats
4428  if (neighbours_txt_file.is_open())
4429  {
4430  neighbours_txt_file<< Direct_string[direction]
4431  << " neighbour of "
4432  << el_pt->object_pt()->number() << " is "
4433  << neighb_pt->object_pt()->number()
4434  << " diff_level " << diff_level
4435  << ". Inside the neighbour the edge is "
4436  << Direct_string[edge] << std::endl;
4437  }
4438 
4439  // Plot neighbour in the appropriate colour
4440  if (neighbours_file.is_open())
4441  {
4442  neighbours_file << "#---------------------------------" << std::endl;
4443  neighbours_file << "#Neighbour element: " << Direct_string[direction]
4444  << std::endl;
4445  neighbours_file << "#---------------------------------" << std::endl;
4446  dynamic_cast<RefineableQElement<3>*>(
4447  neighb_pt->object_pt())->
4448  output_corners(neighbours_file,Colour[direction]);
4449  }
4450 
4451  // Check that local coordinates in the larger element
4452  // lead to the same spatial point as the node vertices
4453  // in the current element
4454  if (neighbours_file.is_open())
4455  {
4456  neighbours_file << "ZONE I=2 C=" << Colour[direction] << std::endl;
4457  }
4458 
4459  // "Low" vertex in the edge
4460  //-------------------------
4461 
4462  // Get coordinates in large (neighbour) element
4463  s[0]=s_lo[0];
4464  s[1]=s_lo[1];
4465  s[2]=s_lo[2];
4466  neighb_pt->object_pt()->get_x(s,x_large);
4467 
4468  // Get coordinates in small element
4469  Vector<double> s(3);
4470  s[0]=S_base_edge(0,direction);
4471  s[1]=S_base_edge(1,direction);
4472  s[2]=S_base_edge(2,direction);
4473  el_pt->object_pt()->get_x(s,x_small);
4474 
4475  double error=
4476  sqrt(pow(x_small[0]-x_large[0],2)+
4477  pow(x_small[1]-x_large[1],2)+
4478  pow(x_small[2]-x_large[2],2));
4479 
4480  if (neighbours_txt_file.is_open())
4481  {
4482  neighbours_txt_file << "Error (1) " << error << std::endl;
4483  }
4484 
4485  if (std::fabs(error)>max_error)
4486  {
4487  max_error=std::fabs(error);
4488  }
4489 
4490 
4491  // Check error and doc mismatch if required
4492  bool stop=false;
4493  std::ofstream mismatch_file;
4494  if (std::fabs(error)>Max_neighbour_finding_tolerance)
4495  {
4496  stop=true;
4497  mismatch_file.open("mismatch.dat");
4498  mismatch_file << "ZONE" << std::endl;
4499  mismatch_file<< x_large[0]<<" "<<x_large[1]<<" "<<x_large[2]<<" 2 \n";
4500  mismatch_file<< x_small[0]<<" "<<x_small[1]<<" "<<x_small[2]<<" 3 \n";
4501  }
4502 
4503  if (neighbours_file.is_open())
4504  {
4505  neighbours_file << "#LOW VERTEX: " << std::endl;
4506  neighbours_file << x_large[0] << " "
4507  << x_large[1] << " "
4508  << x_large[2] <<" 40 \n";
4509  }
4510 
4511 
4512  // "Hi" vertex in the edge
4513  //------------------------
4514 
4515  // Get coordinates in large (neighbour) element
4516  s[0]=s_hi[0];
4517  s[1]=s_hi[1];
4518  s[2]=s_hi[2];
4519  neighb_pt->object_pt()->get_x(s,x_large);
4520 
4521  //Get coordinates in small element
4522  s[0]=S_base_edge(0,direction)+S_step_edge(0,direction);
4523  s[1]=S_base_edge(1,direction)+S_step_edge(1,direction);
4524  s[2]=S_base_edge(2,direction)+S_step_edge(2,direction);
4525  el_pt->object_pt()->get_x(s,x_small);
4526 
4527  // Output
4528  if (neighbours_file.is_open())
4529  {
4530  neighbours_file << "#HI VERTEX: " << std::endl;
4531  neighbours_file << x_large[0]
4532  << " " << x_large[1]
4533  << " " << x_large[2]
4534  << " 80 \n";
4535  }
4536 
4537 
4538  error=
4539  sqrt(pow(x_small[0]-x_large[0],2)+
4540  pow(x_small[1]-x_large[1],2)+
4541  pow(x_small[2]-x_large[2],2));
4542 
4543  if (neighbours_txt_file.is_open())
4544  {
4545  neighbours_txt_file << "Error (2) " << error << std::endl;
4546  }
4547 
4548  if (std::fabs(error)>max_error)
4549  {
4550  max_error=std::fabs(error);
4551  }
4552 
4553 
4554  // Check error and doc mismatch if required
4555  if (std::fabs(error)>Max_neighbour_finding_tolerance) stop=true;
4556  if (stop)
4557  {
4558  if (!mismatch_file.is_open())
4559  {
4560  mismatch_file.open("mismatch.dat");
4561  }
4562  mismatch_file << "ZONE" << std::endl;
4563  mismatch_file << x_large[0] << " " << x_large[1] << " "
4564  << x_large[2] <<" 2 \n";
4565  mismatch_file << x_small[0] << " " << x_small[1] << " "
4566  << x_small[2] <<" 3 \n";
4567  mismatch_file.close();
4568  oomph_info << "Error for direction " << Direct_string[direction]
4569  << " with element "<<i+1<< std::endl;
4570  pause("Error");
4571  }
4572  }
4573 
4574  // If neighbour does not exist: Insert blank zones into file
4575  // so that tecplot can find twelve neighbours for every element
4576  else
4577  {
4578  if (neighbours_file.is_open())
4579  {
4580 
4581  neighbours_file << "#---------------------------------" << std::endl;
4582  neighbours_file << "#no Neighbour in direct: "
4583  << Direct_string[direction]
4584  << std::endl;
4585  neighbours_file << "#---------------------------------" << std::endl;
4586 
4587  dynamic_cast<RefineableQElement<3>*>(
4588  el_pt->object_pt())->output_corners(neighbours_file,"WHITE");
4589  neighbours_file << "ZONE I=2 \n";
4590  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4591  neighbours_file << "-0.05 -0.05 -0.05 0 \n";
4592 
4593 
4594 
4595  // Doc edge for which no neighbour exists but only
4596  // for the smallest elements
4597  if (el_pt->is_leaf())
4598  {
4599  // Get start coordinates in small element
4600  Vector<double> s(3),x_start(3),x_end(3);
4601  s[0]=S_base_edge(0,direction);
4602  s[1]=S_base_edge(1,direction);
4603  s[2]=S_base_edge(2,direction);
4604  el_pt->object_pt()->get_x(s,x_start);
4605 
4606 
4607  //Get coordinates in small element
4608  s[0]=S_base_edge(0,direction)+S_step_edge(0,direction);
4609  s[1]=S_base_edge(1,direction)+S_step_edge(1,direction);
4610  s[2]=S_base_edge(2,direction)+S_step_edge(2,direction);
4611  el_pt->object_pt()->get_x(s,x_end);
4612 
4613  no_true_edge_file << "ZONE I=2" << std::endl;
4614  no_true_edge_file << x_start[0] << " "
4615  << x_start[1] << " "
4616  << x_start[2] << " "
4617  << std::endl;
4618  no_true_edge_file << x_end[0] << " "
4619  << x_end[1] << " "
4620  << x_end[2] << " "
4621  << std::endl;
4622 
4623  }
4624  }
4625 
4626  if (neighbours_file.is_open())
4627  {
4628  neighbours_file << std::endl << std::endl << std::endl;
4629  }
4630  }
4631  } //End of case when element can be documented
4632 
4633  }
4634  }
4635 
4636 }
4637 
4638 
4639 
4640 //////////////////////////////////////////////////////////////////////
4641 //////////////////////////////////////////////////////////////////////
4642 //////////////////////////////////////////////////////////////////////
4643 
4644 
4645 
4646 //================================================================
4647 /// Constructor for OcTreeForest:
4648 ///
4649 /// Pass:
4650 /// - trees_pt[], the Vector of pointers to the constituent trees
4651 /// (OcTreeRoot objects)
4652 ///
4653 //=================================================================
4655  TreeForest(trees_pt)
4656 {
4657 #ifdef LEAK_CHECK
4659 #endif
4660 
4661  // Don't setup neighbours etc. if forest is empty
4662  if (trees_pt.size()==0)
4663  {
4664  return;
4665  }
4666 
4667  using namespace OcTreeNames;
4668 
4669  //Construct the neighbour and rotation scheme, note that all neighbour
4670  //pointers must be set before the constructor is called
4671 
4672 // MemoryUsage::doc_memory_usage(
4673 // "before find_neighbours in octree forest constr");
4674 
4675  // setup the neighbour scheme
4676  find_neighbours();
4677 
4678 // MemoryUsage::doc_memory_usage(
4679 // "after find_neighbours in octree forest constr");
4680 
4681  // setup the rotation scheme
4683 
4684 // MemoryUsage::doc_memory_usage(
4685 // "after construct_up_right_equivalents in octree forest constr");
4686 
4687 }
4688 
4689 //================================================================
4690 /// setup the neighbour scheme : tells to each element in the
4691 /// forest which (if any) element is its {R,L,U,D,B,F} face neighbour
4692 /// and which is its {LB,RB,...,UF} edge neighbour.
4693 //================================================================
4695 {
4696  using namespace OcTreeNames ;
4697  unsigned numtrees=ntree();
4698  int n=0; // to store nnode1d
4699  if(numtrees>0)
4700  {
4701  n=Trees_pt[0]->object_pt()->nnode_1d();
4702  }
4703  else
4704  {
4705  throw OomphLibError(
4706  "Trying to setup the neighbour scheme for an empty forest",
4707  OOMPH_CURRENT_FUNCTION,
4708  OOMPH_EXCEPTION_LOCATION);
4709  }
4710 
4711 
4712  // Number of vertex nodes: 8
4713  unsigned n_vertex_node=8;
4714 
4715  // Find potentially connected trees by identifying
4716  // those whose associated elements share a common vertex node
4717  std::map<Node*,std::set<unsigned> > tree_assoc_with_vertex_node;
4718 
4719  //Loop over all trees
4720  for(unsigned i=0;i<numtrees;i++)
4721  {
4722  // Loop over the vertex nodes of the associated element
4723  for (unsigned j=0;j<n_vertex_node;j++)
4724  {
4725  Node* nod_pt=dynamic_cast<BrickElementBase*>(Trees_pt[i]->object_pt())->
4726  vertex_node_pt(j);
4727  tree_assoc_with_vertex_node[nod_pt].insert(i);
4728  }