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