quadtree.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// Non-inline and non-templated functions for QuadTree and QuadTreeForest
27// classes
28
29// Config header generated by autoconfig
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34#include <set>
35
36// oomph-lib headers
37#include "quadtree.h"
39
40
41namespace oomph
42{
43 //========================================================================
44 /// Bool indicating that static member data has been setup
45 //========================================================================
47
48 //========================================================================
49 /// Colours for neighbours in various directions (static data)
50 //========================================================================
51 Vector<std::string> QuadTree::Colour;
52
53 //=======================================================================
54 /// S_base(i,direction): Initial value for coordinate s[i] on
55 /// the edge indicated by direction (S/E/N/W) (static data)
56 //======================================================================
57 DenseMatrix<double> QuadTree::S_base;
58
59 //======================================================================
60 /// S_step(i,direction) Increments for coordinate s[i] when
61 /// progressing along the edge indicated by direction (S/E/N/W);
62 /// Left/lower vertex: S_base; Right/upper vertex: S_base + S_step
63 /// (static data)
64 //=====================================================================
65 DenseMatrix<double> QuadTree::S_step;
66
67 //=====================================================================
68 /// Translate (enumerated) directions into strings (static data)
69 //=====================================================================
70 Vector<std::string> QuadTree::Direct_string;
71
72 //====================================================================
73 /// Get opposite edge, e.g. Reflect_edge[N]=S (static data)
74 //====================================================================
75 Vector<int> QuadTree::Reflect_edge;
76
77 //===================================================================
78 /// Array of direction/quadrant adjacency scheme:
79 /// Is_adjacent(i_vertex_or_edge,j_quadrant): Is edge/vertex
80 /// adjacent to quadrant? (static data)
81 //===================================================================
82 DenseMatrix<bool> QuadTree::Is_adjacent;
83
84 //====================================================================
85 /// Reflection scheme: Reflect(direction,quadrant): Get mirror
86 /// of quadrant in specified direction. E.g. Reflect(S,NE)=SE
87 /// (static data)
88 //====================================================================
89 DenseMatrix<int> QuadTree::Reflect;
90
91 //=====================================================================
92 /// Rotate coordinates: If north becomes NorthIs then direction
93 /// becomes Rotate(NorthIs,direction). E.g. Rotate(E,NW)=NE;
94 /// (static data)
95 //=====================================================================
96 DenseMatrix<int> QuadTree::Rotate;
97
98 //=====================================================================
99 /// Angle betwen rotated coordinates: If old_direction becomes
100 /// new_direction then the angle between the axes (in anti-clockwise
101 /// direction is Rotate_angle(old_direction,new_direction); E.g.
102 /// Rotate_angle(E,N)=90; (static data)
103 //=====================================================================
104 DenseMatrix<int> QuadTree::Rotate_angle;
105
106 //==========================================================================
107 /// S_direct(direction,son_quadrant): The lower left corner
108 /// of son_quadrant has an offset of h/2 S_direct(direction,son_quadrant)
109 /// in the specified direction. E.g. S_direct(S,NE)=1 and S_direct(S,NW)=0
110 /// (static data)
111 //==========================================================================
112 DenseMatrix<int> QuadTree::S_direct;
113
114
115 //====================================================================
116 /// Setup the static data stored in the QuadTree -- this needs to be
117 /// called before QuadTrees can be used. Automatically called
118 /// by RefineableQuadMesh constructor.
119 //====================================================================
121 {
122 using namespace QuadTreeNames;
123
124
125#ifdef PARANOID
127 {
128 std::ostringstream error_stream;
129 error_stream << "Inconsistent enumeration! \n Tree::OMEGA="
130 << Tree::OMEGA << "\nQuadTree::OMEGA=" << QuadTree::OMEGA
131 << std::endl;
132 throw OomphLibError(
133 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
134 }
135#endif
136
137
138 // Set flag to indicate that static data has been setup
140
141 // Tecplot Colours for neighbours in various directions
142 Colour.resize(27);
143 Colour[SW] = "YELLOW";
144 Colour[SE] = "YELLOW";
145 Colour[NW] = "YELLOW";
146 Colour[NE] = "YELLOW";
147 Colour[E] = "CYAN";
148 Colour[W] = "RED";
149 Colour[N] = "GREEN";
150 Colour[S] = "BLUE";
151 Colour[OMEGA] = "YELLOW";
152
153 // S_base(i,direction): Initial value for coordinate s[i] on the
154 // edge indicated by direction (S/E/N/W)
155 //
156 // S_step(i,direction) Increments for coordinate s[i] when progressing
157 // along the edge indicated by direction (S/E/N/W)
158 S_base.resize(2, 27);
159 S_step.resize(2, 27);
160
161 S_base(0, N) = -1.0;
162 S_base(1, N) = 1.0;
163 S_step(0, N) = 2.0;
164 S_step(1, N) = 0.0;
165
166 S_base(0, S) = -1.0;
167 S_base(1, S) = -1.0;
168 S_step(0, S) = 2.0;
169 S_step(1, S) = 0.0;
170
171 S_base(0, W) = -1.0;
172 S_base(1, W) = -1.0;
173 S_step(0, W) = 0.0;
174 S_step(1, W) = 2.0;
175
176 S_base(0, E) = 1.0;
177 S_base(1, E) = -1.0;
178 S_step(0, E) = 0.0;
179 S_step(1, E) = 2.0;
180
181 // Translate (enumerated) directions into strings
182 Direct_string.resize(27);
183
184 Direct_string[SW] = "SW";
185 Direct_string[NW] = "NW";
186 Direct_string[SE] = "SE";
187 Direct_string[NE] = "NE";
188
189 Direct_string[S] = "S";
190 Direct_string[N] = "N";
191 Direct_string[E] = "E";
192 Direct_string[W] = "W";
193 Direct_string[OMEGA] = "OMEGA";
194
195 // Build direction/quadrant adjacency scheme
196 // Is_adjacent(i_vertex_or_edge,j_quadrant):
197 // Is edge/vertex adjacent to quadrant?
198 Is_adjacent.resize(27, 27);
199
200 Is_adjacent(N, NW) = true;
201 Is_adjacent(E, NW) = false;
202 Is_adjacent(S, NW) = false;
203 Is_adjacent(W, NW) = true;
204 Is_adjacent(NW, NW) = true;
205 Is_adjacent(NE, NW) = false;
206 Is_adjacent(SW, NW) = false;
207 Is_adjacent(SE, NW) = false;
208
209 Is_adjacent(N, NE) = true;
210 Is_adjacent(E, NE) = true;
211 Is_adjacent(S, NE) = false;
212 Is_adjacent(W, NE) = false;
213 Is_adjacent(NW, NE) = false;
214 Is_adjacent(NE, NE) = true;
215 Is_adjacent(SW, NE) = false;
216 Is_adjacent(SE, NE) = false;
217
218 Is_adjacent(N, SW) = false;
219 Is_adjacent(E, SW) = false;
220 Is_adjacent(S, SW) = true;
221 Is_adjacent(W, SW) = true;
222 Is_adjacent(NW, SW) = false;
223 Is_adjacent(NE, SW) = false;
224 Is_adjacent(SW, SW) = true;
225 Is_adjacent(SE, SW) = false;
226
227
228 Is_adjacent(N, SE) = false;
229 Is_adjacent(E, SE) = true;
230 Is_adjacent(S, SE) = true;
231 Is_adjacent(W, SE) = false;
232 Is_adjacent(NW, SE) = false;
233 Is_adjacent(NE, SE) = false;
234 Is_adjacent(SW, SE) = false;
235 Is_adjacent(SE, SE) = true;
236
237 // Rotation scheme: If north becomes NorthIs then
238 // direction becomes Rotate(NorthIs,direction)
239 // Initialise to OMEGA
240 Rotate.resize(27, 27);
241
242 Rotate(N, N) = N;
243 Rotate(N, E) = E;
244 Rotate(N, S) = S;
245 Rotate(N, W) = W;
246 Rotate(N, NW) = NW;
247 Rotate(N, NE) = NE;
248 Rotate(N, SE) = SE;
249 Rotate(N, SW) = SW;
250
251 Rotate(W, N) = W;
252 Rotate(W, E) = N;
253 Rotate(W, S) = E;
254 Rotate(W, W) = S;
255 Rotate(W, NW) = SW;
256 Rotate(W, NE) = NW;
257 Rotate(W, SE) = NE;
258 Rotate(W, SW) = SE;
259
260 Rotate(S, N) = S;
261 Rotate(S, E) = W;
262 Rotate(S, S) = N;
263 Rotate(S, W) = E;
264 Rotate(S, NW) = SE;
265 Rotate(S, NE) = SW;
266 Rotate(S, SE) = NW;
267 Rotate(S, SW) = NE;
268
269 Rotate(E, N) = E;
270 Rotate(E, E) = S;
271 Rotate(E, S) = W;
272 Rotate(E, W) = N;
273 Rotate(E, NW) = NE;
274 Rotate(E, NE) = SE;
275 Rotate(E, SE) = SW;
276 Rotate(E, SW) = NW;
277
278 // Angle betwen rotated coordinates:
279 // old_direction becomes new_direction then the angle between
280 // the axes is Rotate_angle(old_direction,new_direction)
281 Rotate_angle.resize(27, 27);
282
283 Rotate_angle(N, N) = 0;
284 Rotate_angle(N, W) = 90;
285 Rotate_angle(N, S) = 180;
286 Rotate_angle(N, E) = 270;
287
288 Rotate_angle(S, S) = 0;
289 Rotate_angle(S, E) = 90;
290 Rotate_angle(S, N) = 180;
291 Rotate_angle(S, W) = 270;
292
293 Rotate_angle(W, W) = 0;
294 Rotate_angle(W, S) = 90;
295 Rotate_angle(W, E) = 180;
296 Rotate_angle(W, N) = 270;
297
298 Rotate_angle(E, E) = 0;
299 Rotate_angle(E, N) = 90;
300 Rotate_angle(E, W) = 180;
301 Rotate_angle(E, S) = 270;
302
303 // Reflection scheme:
304 // Reflect(direction,quadrant): Get mirror of quadrant in direction
305 Reflect.resize(27, 27);
306
307 Reflect(N, NW) = SW;
308 Reflect(E, NW) = NE;
309 Reflect(S, NW) = SW;
310 Reflect(W, NW) = NE;
311 Reflect(NW, NW) = SE;
312 Reflect(NE, NW) = SE;
313 Reflect(SW, NW) = SE;
314 Reflect(SE, NW) = SE;
315
316 Reflect(N, NE) = SE;
317 Reflect(E, NE) = NW;
318 Reflect(S, NE) = SE;
319 Reflect(W, NE) = NW;
320 Reflect(NW, NE) = SW;
321 Reflect(NE, NE) = SW;
322 Reflect(SW, NE) = SW;
323 Reflect(SE, NE) = SW;
324
325 Reflect(N, SW) = NW;
326 Reflect(E, SW) = SE;
327 Reflect(S, SW) = NW;
328 Reflect(W, SW) = SE;
329 Reflect(NW, SW) = NE;
330 Reflect(NE, SW) = NE;
331 Reflect(SW, SW) = NE;
332 Reflect(SE, SW) = NE;
333
334 Reflect(N, SE) = NE;
335 Reflect(E, SE) = SW;
336 Reflect(S, SE) = NE;
337 Reflect(W, SE) = SW;
338 Reflect(NW, SE) = NW;
339 Reflect(NE, SE) = NW;
340 Reflect(SW, SE) = NW;
341 Reflect(SE, SE) = NW;
342
343 // S_direct(direction,son_quadrant): The lower left corner of
344 // my son_quadrant has an offset of h/2 S_direct(direction,son_quadrant)
345 // in the direction. [Look in the positive direction along the
346 // edge 'direction' and measure the offset of the lower left corner
347 // of the son quadrant in this direction]
348 S_direct.resize(27, 27);
349
350 S_direct(W, NW) = 1;
351 S_direct(E, NW) = 1;
352 S_direct(N, NW) = 0;
353 S_direct(S, NW) = 0;
354
355 S_direct(W, NE) = 1;
356 S_direct(E, NE) = 1;
357 S_direct(N, NE) = 1;
358 S_direct(S, NE) = 1;
359
360 S_direct(W, SW) = 0;
361 S_direct(E, SW) = 0;
362 S_direct(N, SW) = 0;
363 S_direct(S, SW) = 0;
364
365 S_direct(W, SE) = 0;
366 S_direct(E, SE) = 0;
367 S_direct(N, SE) = 1;
368 S_direct(S, SE) = 1;
369
370 // Get opposite edge, e.g. Reflect_edge(N)=S
371 Reflect_edge.resize(27);
372
373 Reflect_edge[N] = S;
374 Reflect_edge[S] = N;
375 Reflect_edge[E] = W;
376 Reflect_edge[W] = E;
377 }
378
379 //================================================================
380 /// Return pointer to greater or equal-sized edge neighbour
381 /// in specified \c direction; also provide info regarding the relative
382 /// size and orientation of neighbour:
383 /// - The vector translate_s turns the index of the local coordinate
384 /// in the present quadtree into that of the neighbour. If there are no
385 /// rotations then translate_s[i] = i, but if, for example, the neighbour's
386 /// eastern face is adjacent to our northern face translate_s[0] = 1
387 /// and translate_s[1] = 0. Of course, this could be deduced after the
388 /// fact, but it's easier to do it here.
389 /// - Let's have a look at the current quadtree's edge in the specified
390 /// direction. The edge will be parallel to one of the two
391 /// local coordinates in the element. This edge has two
392 /// vertices. The vertex at the minimum value of the local
393 /// coordinate (the "lo" vertex) is located
394 /// at the local coordinates (\c s_lo[0], \c s_lo[1]) in the neighbouring
395 /// quadtree.
396 /// - ditto with s_hi: The vertex at the maximum value of the local
397 /// coordinate located at the local coordinates (\c s_hi[0],
398 /// \c s_hi[1]) in the neighbouring quadtree.
399 /// - We're looking for a neighbour in the specified \c direction. When
400 /// viewed from the neighbouring quadtree, the edge that separates
401 /// the present quadtree from its neighbour is the neighbour's \c edge
402 /// edge. If there's no rotation between the two quadtrees, this is a
403 /// simple reflection: For instance, if we're looking
404 /// for a neighhbour in the \c N [orthern] \c direction, \c edge will
405 /// be \c S [outh]
406 /// - \c diff_level <= 0 indicates the difference in refinement levels between
407 /// the two neighbours. If \c diff_level==0, the neighbour has the
408 /// same size as the current quadtree.
409 /// - \c in_neighbouring_tree returns true is we have had to flip
410 /// to a different root, even if that root is actually the same
411 /// as it can be in periodic problems.
412 //=================================================================
414 Vector<unsigned>& translate_s,
415 Vector<double>& s_lo,
416 Vector<double>& s_hi,
417 int& edge,
418 int& diff_level,
419 bool& in_neighbouring_tree) const
420 {
421 using namespace QuadTreeNames;
422
423#ifdef PARANOID
424 if ((direction != S) && (direction != E) && (direction != N) &&
425 (direction != W))
426 {
427 std::ostringstream error_stream;
428 error_stream << "Direction " << direction << " is not N, S, E, W"
429 << std::endl;
430
431 throw OomphLibError(
432 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
433 }
434#endif
435
436 // Initialise in_neighbouring tree to false. It will be set to true
437 // during the recursion if we do actually hop over in to the neighbour
438 in_neighbouring_tree = false;
439
440 // Maximum level to which we're allowed to descend (we only want
441 // greater-or-equal-sized neighbours)
442 int max_level = Level;
443
444 // Current element has the following root:
445 QuadTreeRoot* orig_root_pt = dynamic_cast<QuadTreeRoot*>(Root_pt);
446
447 // Initialise offset in local coordinate
448 double s_diff = 0;
449
450 // Initialise difference in level
451 diff_level = 0;
452
453 // Find neighbour
454 QuadTree* return_pt = gteq_edge_neighbour(direction,
455 s_diff,
456 diff_level,
457 in_neighbouring_tree,
458 max_level,
459 orig_root_pt);
460
461 QuadTree* neighb_pt = return_pt;
462
463 // If neighbour exists: What's the direction of the interfacial
464 // edge when viewed from within the neighbour element?
465 if (neighb_pt != 0)
466 {
467 // Rotate things around (the orientation of N in the neighbour might be
468 // be different from that in the present element)
469 // Initialise the direction
470 int new_dir = direction;
471 // If the neighbour has a different root, then there could be a possible
472 // rotation, find it
473 if (neighb_pt->root_pt() != Root_pt)
474 {
475 new_dir = Rotate(orig_root_pt->north_equivalent(direction), direction);
476 }
477
478 s_lo[0] = S_base(0, Reflect_edge[new_dir]) +
479 S_step(0, Reflect_edge[new_dir]) * s_diff;
480 s_lo[1] = S_base(1, Reflect_edge[new_dir]) +
481 S_step(1, Reflect_edge[new_dir]) * s_diff;
482
483 s_hi[0] = S_base(0, Reflect_edge[new_dir]) +
484 S_step(0, Reflect_edge[new_dir]) * pow(2.0, diff_level) +
485 S_step(0, Reflect_edge[new_dir]) * s_diff;
486 s_hi[1] = S_base(1, Reflect_edge[new_dir]) +
487 S_step(1, Reflect_edge[new_dir]) * pow(2.0, diff_level) +
488 S_step(1, Reflect_edge[new_dir]) * s_diff;
489
490 Vector<double> s_lo_new(2), s_hi_new(2);
491
492 // What's the direction of the interfacial edge when viewed from within
493 // the neighbour element?
494 edge = Reflect_edge[new_dir];
495
496 // Set up the translation scheme for the local coordinates
497 {
498 bool swap = false;
499 // Do we need to switch the coordinate
500 switch (direction)
501 {
502 // If the direction is north or south,
503 // but the neighbour's coordinate s[1] is not constant, we must swap
504 case N:
505 case S:
506 if (s_lo[1] != s_hi[1])
507 {
508 swap = true;
509 }
510 break;
511 // If the direction is east or west
512 // but the neighbour's corodinate s[0] is not constant, we must swap
513 case E:
514 case W:
515 if (s_lo[0] != s_hi[0])
516 {
517 swap = true;
518 }
519 break;
520 // Catch all totally un-necessary
521 default:
522 std::ostringstream error_stream;
523 error_stream << "Direction " << direction << " is not N, S, E, W"
524 << std::endl;
525
526 throw OomphLibError(error_stream.str(),
527 OOMPH_CURRENT_FUNCTION,
528 OOMPH_EXCEPTION_LOCATION);
529 }
530
531 // If we must swap, then do so
532 if (swap)
533 {
534 translate_s[0] = 1;
535 translate_s[1] = 0;
536 }
537 // Othewise, it's just a straight translation
538 else
539 {
540 translate_s[0] = 0;
541 translate_s[1] = 1;
542 }
543 }
544
545 // Reverse directions?
546 s_lo_new[0] = s_lo[0];
547 s_lo_new[1] = s_lo[1];
548 s_hi_new[0] = s_hi[0];
549 s_hi_new[1] = s_hi[1];
550
551 if (((edge == N) || (edge == S)) &&
552 ((Rotate_angle(direction, new_dir) == 90) ||
553 (Rotate_angle(direction, new_dir) == 180)))
554 {
555 s_lo_new[0] = -s_lo[0];
556 s_hi_new[0] = -s_hi[0];
557 }
558 if (((edge == E) || (edge == W)) &&
559 ((Rotate_angle(direction, new_dir) == 270) ||
560 (Rotate_angle(direction, new_dir) == 180)))
561 {
562 s_lo_new[1] = -s_lo[1];
563 s_hi_new[1] = -s_hi[1];
564 }
565
566 s_lo[0] = s_lo_new[0];
567 s_lo[1] = s_lo_new[1];
568 s_hi[0] = s_hi_new[0];
569 s_hi[1] = s_hi_new[1];
570 }
571 return return_pt;
572 }
573
574 //================================================================
575 /// Find `greater-or-equal-sized edge neighbour' in given direction
576 /// (N/E/S/W).
577 ///
578 /// This is an auxiliary routine which allows neighbour finding in adjacent
579 /// quadtrees. Needs to keep track of previous son types and
580 /// the maximum level to which search is performed.
581 ///
582 /// Parameters:
583 ///
584 /// - direction: N/S/E/W: Direction in which neighbour has to be found.
585 /// - s_diff: Offset of lower/left vertex from corresponding vertex in
586 /// neighbour. Note that this is input/output as it needs to be incremented/
587 /// decremented during the recursive calls to this function.
588 /// - edge: We're looking for the neighbour across our edge 'direction'
589 /// (N/S/E/W). When viewed from the neighbour, this edge is
590 /// `edge' (N/S/E/W). [If there's no relative rotation between neighbours
591 /// then this is a mere reflection, e.g. direction=N --> edge=S etc.]
592 /// - diff_level <= 0 indicates the difference in quadtree levels
593 /// between the current element and its neighbour.
594 /// - max_level is the maximum level to which the neighbour search is
595 /// allowed to proceed. This is again necessary because in a forest,
596 /// the neighbour search isn't based on pure recursion.
597 /// - orig_root_pt identifies the root node of the element whose
598 /// neighbour we're really trying to find by all these recursive calls.
599 ///
600 //=================================================================
602 const int& direction,
603 double& s_diff,
604 int& diff_level,
605 bool& in_neighbouring_tree,
606 int max_level,
607 QuadTreeRoot* const& orig_root_pt) const
608 {
609 using namespace QuadTreeNames;
610
611#ifdef PARANOID
612 if ((direction != S) && (direction != E) && (direction != N) &&
613 (direction != W))
614 {
615 std::ostringstream error_stream;
616 error_stream << "Direction " << direction << " is not N, S, E, W"
617 << std::endl;
618
619 throw OomphLibError(
620 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
621 }
622#endif
623
624 QuadTree* next_el_pt;
625 QuadTree* return_el_pt;
626
627 // STEP 1: Find the neighbour's father
628 //--------
629
630 // Does the element have a father?
631 if (Father_pt != 0)
632 {
633 // If the present quadrant (whose location inside its
634 // father element is specified by Son_type) is adjacent to the
635 // father's edge in the required direction, then its neighbour has
636 // a different father ---> we need to climb up the tree to
637 // the father and find his neighbour in the required direction
638
639 // Note that this is the cunning recursive part. The returning may not
640 // stop until we hit the very top of the tree, when the element does NOT
641 // have a father
642 if (Is_adjacent(direction, Son_type))
643 {
644 next_el_pt = dynamic_cast<QuadTree*>(Father_pt)->gteq_edge_neighbour(
645 direction,
646 s_diff,
647 diff_level,
648 in_neighbouring_tree,
649 max_level,
650 orig_root_pt);
651 }
652 // If the present quadrant is not adjacent to the
653 // father's edge in the required direction, then the
654 // neighbour has the same father and is obtained
655 // by the appropriate reflection inside the father element
656 // This will only be called if we have not left the original tree.
657 else
658 {
659 next_el_pt = dynamic_cast<QuadTree*>(Father_pt);
660 }
661
662 // We're about to ascend one level:
663 diff_level -= 1;
664
665 // Work out position of lower (or left) corner of present edge
666 // in its father element
667 s_diff += pow(0.5, -diff_level) * S_direct(direction, Son_type);
668
669 // STEP 2: We have now located the neighbour's father and need to
670 // ------- find the appropriate son.
671
672 // Buffer cases where the neighbour (and hence its father) lie outside
673 // the boundary
674 if (next_el_pt != 0)
675 {
676 // If the father is a leaf then we can't descend to the same
677 // level as the present node ---> simply return the father himself
678 // as the (greater) neighbour. Same applies if we are about
679 // to descend lower than the max_level (in a neighbouring tree)
680 if ((next_el_pt->Son_pt.size() == 0) ||
681 (next_el_pt->Level > max_level - 1))
682 {
683 return_el_pt = next_el_pt;
684 }
685 // We have located the neighbour's father: The position of the
686 // neighbour is obtained by `reflecting' the position of the
687 // node itself.
688
689 // We know exactly how to reflect, because we know which son type we
690 // are and we have the pointer to the neighbours father
691 else
692 {
693 int son_quadrant = Reflect(direction, Son_type);
694
695 // If the root of the neighbour's father is not our root, we
696 // might need to rotate
697 if (orig_root_pt != next_el_pt->Root_pt)
698 {
699 // Get the north equivalent of the next element
700 int my_north =
701 dynamic_cast<QuadTreeRoot*>(Root_pt)->north_equivalent(direction);
702 son_quadrant = Rotate(my_north, son_quadrant);
703 }
704
705 // The next element in the tree is the appropriate son of the
706 // neighbour's father
707 return_el_pt =
708 dynamic_cast<QuadTree*>(next_el_pt->Son_pt[son_quadrant]);
709
710 // Work out position of lower (or left) corner of present edge
711 // in next higher element
712 s_diff -= pow(0.5, -diff_level) * S_direct(direction, Son_type);
713
714 // We have just descended one level
715 diff_level += 1;
716 }
717 }
718 // The neighbour's father lies outside the boundary --> the neighbour
719 // itself does too --> return NULL.
720 else
721 {
722 return_el_pt = 0;
723 }
724 }
725 // Element does not have a father --> check if it has a neighbouring
726 // tree in the appropriate direction
727 else
728 {
729 // Find neighbouring root
730 if (Root_pt->neighbour_pt(direction) != 0)
731 {
732 // In this case we have moved to a neighbour, so set the flag
733 in_neighbouring_tree = true;
734 return_el_pt =
735 dynamic_cast<QuadTreeRoot*>(Root_pt->neighbour_pt(direction));
736 }
737 // No neighbouring tree, so there really is no neighbour --> return NULL
738 else
739 {
740 return_el_pt = 0;
741 }
742 }
743
744 return return_el_pt;
745 }
746
747 //================================================================
748 /// Traverse Tree: Preorder traverse and stick pointers to
749 /// neighbouring leaf nodes (only) into Vector
750 //=================================================================
752 Vector<const QuadTree*>& tree_neighbouring_nodes,
753 Vector<Vector<double>>& tree_neighbouring_s_lo,
754 Vector<Vector<double>>& tree_neighbouring_s_hi,
755 Vector<int>& tree_neighbouring_diff_level,
756 const QuadTree* my_neigh_pt,
757 const int& direction) const
758 {
759 // If the tree has sons
760 unsigned numsons = Son_pt.size();
761 if (numsons > 0)
762 {
763 // Now do the sons (if they exist)
764 for (unsigned i = 0; i < numsons; i++)
765 {
766 dynamic_cast<QuadTree*>(Son_pt[i])
767 ->stick_neighbouring_leaves_into_vector(tree_neighbouring_nodes,
768 tree_neighbouring_s_lo,
769 tree_neighbouring_s_hi,
770 tree_neighbouring_diff_level,
771 my_neigh_pt,
772 direction);
773 }
774 }
775 else
776 {
777 // Required data for neighbour-finding routine
778 Vector<unsigned> translate_s(2);
779 Vector<double> s_lo(2), s_hi(2);
780 int edge, diff_level;
781 bool in_neighbouring_tree;
782 QuadTree* neigh_pt;
783
784 // Get neighbouring tree
785 neigh_pt = gteq_edge_neighbour(direction,
786 translate_s,
787 s_lo,
788 s_hi,
789 edge,
790 diff_level,
791 in_neighbouring_tree);
792
793 // Check if the neighbour is the same as the tree passed in
794 // (i.e. Am I a neighbour of the master element's tree?)
795 if (neigh_pt == my_neigh_pt)
796 {
797 // Add the element and the diff_level to passed vectors
798 tree_neighbouring_nodes.push_back(this);
799 tree_neighbouring_s_lo.push_back(s_lo);
800 tree_neighbouring_s_hi.push_back(s_hi);
801 tree_neighbouring_diff_level.push_back(diff_level);
802 }
803 }
804 }
805
806 //================================================================
807 /// Self-test: Check neighbour finding routine. For each element
808 /// in the tree and for each vertex, determine the
809 /// distance between the vertex and its position in the
810 /// neigbour. If the difference is less than
811 /// Tree::Max_neighbour_finding_tolerance.
812 /// return success (0), otherwise failure (1)
813 //=================================================================
815 {
816 // Stick pointers to all nodes into Vector and number elements
817 // in the process
818 Vector<Tree*> all_nodes_pt;
820 long int count = 0;
821 unsigned long num_nodes = all_nodes_pt.size();
822 for (unsigned long i = 0; i < num_nodes; i++)
823 {
824 all_nodes_pt[i]->object_pt()->set_number(++count);
825 }
826
827 // Check neighbours (distance between hanging nodes) -- don't print (keep
828 // output streams closed)
829 double max_error = 0.0;
830 std::ofstream neighbours_file;
831 std::ofstream neighbours_txt_file;
833 all_nodes_pt, neighbours_file, neighbours_txt_file, max_error);
834
835 if (max_error > Max_neighbour_finding_tolerance)
836 {
837 oomph_info << "\n \n Failed self_test() for QuadTree: Max. error "
838 << max_error << std::endl
839 << std::endl;
840 return 1;
841 }
842 else
843 {
844 oomph_info << "\n \n Passed self_test() for QuadTree: Max. error "
845 << max_error << std::endl
846 << std::endl;
847 return 0;
848 }
849 }
850
851
852 //================================================================
853 /// Constructor for QuadTreeForest:
854 ///
855 /// Pass:
856 /// - trees_pt[], the Vector of pointers to the constituent trees
857 /// (QuadTreeRoot objects)
858 ///
859 /// Note that the pointers to the neighbour's of each tree must have
860 /// been allocated before the constructor is called, otherwise the
861 /// relative rotation scheme will not be constructed correctly.
862 //=================================================================
864 : TreeForest(trees_pt)
865 {
866#ifdef LEAK_CHECK
868#endif
869
870 // Don't setup neighbours etc. if forest is empty
871 if (trees_pt.size() == 0)
872 {
873 return;
874 }
875
876 using namespace QuadTreeNames;
877
878 // Setup the neighbours
880
881 // Construct the rotation scheme, note that all neighbour pointers must
882 // be set before the constructor is called
884 }
885
886
887 //================================================================
888 /// Setup the neighbour lookup schemes for all constituent
889 /// quadtrees.
890 //================================================================
892 {
893 using namespace QuadTreeNames;
894
895 unsigned numtrees = ntree();
896 unsigned n = 0; // to store nnode1d
897 if (numtrees > 0)
898 {
899 n = Trees_pt[0]->object_pt()->nnode_1d();
900 }
901 else
902 {
903 throw OomphLibError(
904 "Trying to setup the neighbour scheme for an empty forest\n",
905 OOMPH_CURRENT_FUNCTION,
906 OOMPH_EXCEPTION_LOCATION);
907 }
908
909 // Number of vertex nodes: 4
910 unsigned n_vertex_node = 4;
911
912 // Find potentially connected trees by identifying
913 // those whose associated elements share a common vertex node
914 std::map<Node*, std::set<unsigned>> tree_assoc_with_vertex_node;
915
916 // Loop over all trees
917 for (unsigned i = 0; i < numtrees; i++)
918 {
919 // Loop over the vertex nodes of the associated element
920 for (unsigned j = 0; j < n_vertex_node; j++)
921 {
922 Node* nod_pt = dynamic_cast<QuadElementBase*>(Trees_pt[i]->object_pt())
923 ->vertex_node_pt(j);
924 tree_assoc_with_vertex_node[nod_pt].insert(i);
925 }
926 }
927
928
929 // For each tree we store a set of potentially neighbouring trees
930 // i.e. trees that share at least one node
931 Vector<std::set<unsigned>> potentially_neighb_tree(numtrees);
932
933 // Loop over vertex nodes
934 for (std::map<Node*, std::set<unsigned>>::iterator it =
935 tree_assoc_with_vertex_node.begin();
936 it != tree_assoc_with_vertex_node.end();
937 it++)
938 {
939 // Loop over connected elements twice
940 for (std::set<unsigned>::iterator it_el1 = it->second.begin();
941 it_el1 != it->second.end();
942 it_el1++)
943 {
944 unsigned i = (*it_el1);
945 for (std::set<unsigned>::iterator it_el2 = it->second.begin();
946 it_el2 != it->second.end();
947 it_el2++)
948 {
949 unsigned j = (*it_el2);
950 // These two elements are potentially connected
951 if (i != j)
952 {
953 potentially_neighb_tree[i].insert(j);
954 }
955 }
956 }
957 }
958
959
960 // Loop over all trees
961 for (unsigned i = 0; i < numtrees; i++)
962 {
963 // Loop over their potential neighbours
964 for (std::set<unsigned>::iterator it = potentially_neighb_tree[i].begin();
965 it != potentially_neighb_tree[i].end();
966 it++)
967 {
968 unsigned j = (*it);
969
970 // is it the Northern neighbour ?
971 bool is_N_neighbour =
972 ((Trees_pt[j]->object_pt()->get_node_number(
973 Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1) &&
974 (Trees_pt[j]->object_pt()->get_node_number(
975 Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
976
977
978 // is it the Southern neighbour ?
979 bool is_S_neighbour =
980 ((Trees_pt[j]->object_pt()->get_node_number(
981 Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
982 (Trees_pt[j]->object_pt()->get_node_number(
983 Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1));
984
985
986 // is it the Eastern neighbour ?
987 bool is_E_neighbour =
988 ((Trees_pt[j]->object_pt()->get_node_number(
989 Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1) &&
990 (Trees_pt[j]->object_pt()->get_node_number(
991 Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
992
993 // is it the Western neighbour ?
994 bool is_W_neighbour =
995 ((Trees_pt[j]->object_pt()->get_node_number(
996 Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
997 (Trees_pt[j]->object_pt()->get_node_number(
998 Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1));
999
1000
1001 if (is_N_neighbour) Trees_pt[i]->neighbour_pt(N) = Trees_pt[j];
1002 if (is_S_neighbour) Trees_pt[i]->neighbour_pt(S) = Trees_pt[j];
1003 if (is_E_neighbour) Trees_pt[i]->neighbour_pt(E) = Trees_pt[j];
1004 if (is_W_neighbour) Trees_pt[i]->neighbour_pt(W) = Trees_pt[j];
1005 }
1006 }
1007
1008
1009 // Old hacky version with horrendous scaling
1010 if (false)
1011 {
1012 // Loop over all trees
1013 for (unsigned i = 0; i < numtrees; i++)
1014 {
1015 // Loop over all the other elements
1016 for (unsigned j = 0; j < numtrees; j++)
1017 {
1018 if (j != i) // make sure we are not looking at the element itself
1019 {
1020 // is it the Northern neighbour ?
1021 bool is_N_neighbour =
1022 ((Trees_pt[j]->object_pt()->get_node_number(
1023 Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1) &&
1024 (Trees_pt[j]->object_pt()->get_node_number(
1025 Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
1026
1027
1028 // is it the Southern neighbour ?
1029 bool is_S_neighbour =
1030 ((Trees_pt[j]->object_pt()->get_node_number(
1031 Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
1032 (Trees_pt[j]->object_pt()->get_node_number(
1033 Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1));
1034
1035
1036 // is it the Eastern neighbour ?
1037 bool is_E_neighbour =
1038 ((Trees_pt[j]->object_pt()->get_node_number(
1039 Trees_pt[i]->object_pt()->node_pt(n - 1)) != -1) &&
1040 (Trees_pt[j]->object_pt()->get_node_number(
1041 Trees_pt[i]->object_pt()->node_pt(n * n - 1)) != -1));
1042
1043 // is it the Western neighbour ?
1044 bool is_W_neighbour =
1045 ((Trees_pt[j]->object_pt()->get_node_number(
1046 Trees_pt[i]->object_pt()->node_pt(0)) != -1) &&
1047 (Trees_pt[j]->object_pt()->get_node_number(
1048 Trees_pt[i]->object_pt()->node_pt(n * (n - 1))) != -1));
1049
1050
1051 if (is_N_neighbour) Trees_pt[i]->neighbour_pt(N) = Trees_pt[j];
1052 if (is_S_neighbour) Trees_pt[i]->neighbour_pt(S) = Trees_pt[j];
1053 if (is_E_neighbour) Trees_pt[i]->neighbour_pt(E) = Trees_pt[j];
1054 if (is_W_neighbour) Trees_pt[i]->neighbour_pt(W) = Trees_pt[j];
1055 }
1056 }
1057 }
1058 }
1059 }
1060
1061
1062 //================================================================
1063 /// Construct the rotation scheme for the quadtree forest.
1064 /// Note that all pointers to neighbours must have been allocated
1065 /// for this to work.
1066 //================================================================
1068 {
1069 using namespace QuadTreeNames;
1070
1071 unsigned numtrees = ntree();
1072 // Loop over all the trees
1073 for (unsigned i = 0; i < numtrees; i++)
1074 {
1075 // Find the pointer to the northern neighbour
1076 QuadTreeRoot* neigh_pt = quad_neigh_pt(i, N);
1077 // If there is a neighbour
1078 if (neigh_pt != 0)
1079 {
1080 // Find the direction of the present tree, as viewed from the neighbour
1081 int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
1082
1083 // Set up the rotation scheme
1084 switch (direction)
1085 {
1086 // If N neighbour has this tree on S, north equivalent is N
1087 case S:
1089 break;
1090 // If N neighbour has this tree on W, north equivalent is E
1091 case W:
1093 break;
1094 // If N neighbour has this tree on N, north equivalent is S
1095 case N:
1097 break;
1098 // If N neighbour has this tree on E, north equivalent is W
1099 case E:
1101 break;
1102 // If N neighbour does not have pointer to this tree, die
1103 default:
1104 std::ostringstream error_stream;
1105 error_stream
1106 << "Tree " << i
1107 << "'s Northern neighbour has no neighbour pointer to Tree " << i
1108 << std::endl;
1109
1110 throw OomphLibError(error_stream.str(),
1111 OOMPH_CURRENT_FUNCTION,
1112 OOMPH_EXCEPTION_LOCATION);
1113 }
1114 }
1115
1116 // Find the pointer to the eastern neighbour
1117 neigh_pt = quad_neigh_pt(i, E);
1118 // If there is a neighbour
1119 if (neigh_pt != 0)
1120 {
1121 // Find the direction of the present tree, as viewed from the neighbour
1122 int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
1123
1124 // Set up the rotation scheme
1125 switch (direction)
1126 {
1127 // If E neighbour has this tree on W, north equivalent is N
1128 case W:
1130 break;
1131 // If E neighbour has this tree on N, north equivalent is E
1132 case N:
1134 break;
1135 // If E neighbour has this tree on E, north equivalent is S
1136 case E:
1138 break;
1139 // If E neighbour has this tree on S, north equivalent is W
1140 case S:
1142 break;
1143 // If E neighbour does not have pointer to this tree, die
1144 default:
1145 std::ostringstream error_stream;
1146 error_stream
1147 << "Tree " << i
1148 << "'s Eastern neighbour has no neighbour pointer to Tree " << i
1149 << std::endl;
1150
1151 throw OomphLibError(error_stream.str(),
1152 OOMPH_CURRENT_FUNCTION,
1153 OOMPH_EXCEPTION_LOCATION);
1154 }
1155 }
1156
1157 // Find the pointer to the southern neighbour
1158 neigh_pt = quad_neigh_pt(i, S);
1159 // If there is a neighbour
1160 if (neigh_pt != 0)
1161 {
1162 // Find the direction of the present tree, as viewed from the neighbour
1163 int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
1164
1165 // Set up the rotation scheme
1166 switch (direction)
1167 {
1168 // If S neighbour has this tree on N, north equivalent is N
1169 case N:
1171 break;
1172 // If S neighbour has this tree on E, north equivalent is E
1173 case E:
1175 break;
1176 // If S neighbour has this tree on S, north equivalent is S
1177 case S:
1179 break;
1180 // If S neighbour has this tree on W, north equivalent is W
1181 case W:
1183 break;
1184 // If S neighbour does not have pointer to this tree, die
1185 default:
1186 std::ostringstream error_stream;
1187 error_stream
1188 << "Tree " << i
1189 << "'s Southern neighbour has no neighbour pointer to Tree " << i
1190 << std::endl;
1191
1192 throw OomphLibError(error_stream.str(),
1193 OOMPH_CURRENT_FUNCTION,
1194 OOMPH_EXCEPTION_LOCATION);
1195 }
1196 }
1197
1198 // Find the pointer to the western neighbour
1199 neigh_pt = quad_neigh_pt(i, W);
1200 // If there is a neighbour
1201 if (neigh_pt != 0)
1202 {
1203 // Find the direction of the present tree, as viewed from the neighbour
1204 int direction = neigh_pt->direction_of_neighbour(quadtree_pt(i));
1205
1206 // Set up the rotation scheme
1207 switch (direction)
1208 {
1209 // If W neighbour has this tree on E, north equivalent is N
1210 case E:
1212 break;
1213 // If W neighbour has this tree on S, north equivalent is E
1214 case S:
1216 break;
1217 // If W neighbour has this tree on W, north equivalent is S
1218 case W:
1220 break;
1221 // If W neighbour has this tree on N, north equivalent is W
1222 case N:
1224 break;
1225 // If W neighbour does not have pointer to this tree, die
1226 default:
1227 std::ostringstream error_stream;
1228 error_stream
1229 << "Tree " << i
1230 << "'s Western neighbour has no neighbour pointer to Tree " << i
1231 << std::endl;
1232
1233 throw OomphLibError(error_stream.str(),
1234 OOMPH_CURRENT_FUNCTION,
1235 OOMPH_EXCEPTION_LOCATION);
1236 }
1237 }
1238 }
1239 }
1240
1241 //================================================================
1242 /// Document and check all the neighbours in all the nodes in
1243 /// the forest
1244 //================================================================
1246 {
1247 // Create Vector of elements
1248 Vector<Tree*> all_tree_nodes_pt;
1249 this->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
1250
1251 // Create storage for information files
1252 std::ofstream neigh_file;
1253 std::ofstream neigh_txt_file;
1254
1255 // If we are documenting the results, then open the files
1256 if (doc_info.is_doc_enabled())
1257 {
1258 std::ostringstream fullname;
1259 fullname << doc_info.directory() << "/neighbours" << doc_info.number()
1260 << ".dat";
1261 oomph_info << "opened " << fullname.str() << " to doc neighbours"
1262 << std::endl;
1263 neigh_file.open(fullname.str().c_str());
1264 fullname.str("");
1265 fullname << doc_info.directory() << "/neighbours" << doc_info.number()
1266 << ".txt";
1267 oomph_info << "opened " << fullname.str() << " to doc neighbours"
1268 << std::endl;
1269 neigh_txt_file.open(fullname.str().c_str());
1270 }
1271
1272 // Call the standard documentation function
1273 double max_error = 0.0;
1275 all_tree_nodes_pt, neigh_file, neigh_txt_file, max_error);
1276
1277 // If the error is too large complain
1278 if (max_error > Tree::max_neighbour_finding_tolerance())
1279 {
1280 std::ostringstream error_stream;
1281 error_stream << "Max. error in quadtree neighbour finding: " << max_error
1282 << " is too big" << std::endl;
1283 error_stream
1284 << "i.e. bigger than Tree::max_neighbour_finding_tolerance()="
1285 << Tree::max_neighbour_finding_tolerance() << std::endl;
1286
1287 // Close the files if they were opened
1288 if (doc_info.is_doc_enabled())
1289 {
1290 neigh_file.close();
1291 neigh_txt_file.close();
1292 }
1293
1294 throw OomphLibError(
1295 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1296 }
1297 else
1298 {
1299 oomph_info << "Max. error in quadtree neighbour finding: " << max_error
1300 << " is OK" << std::endl;
1302 << "i.e. less than QuadTree::max_neighbour_finding_tolerance()="
1304 }
1305
1306 // Close the files if they were opened
1307 if (doc_info.is_doc_enabled())
1308 {
1309 neigh_file.close();
1310 neigh_txt_file.close();
1311 }
1312 }
1313
1314 //================================================================
1315 /// Open output files that will stored any hanging nodes that are
1316 /// created in the mesh refinement process.
1317 /// ===============================================================
1319 DocInfo& doc_info, Vector<std::ofstream*>& output_stream)
1320 {
1321 // In 2D, there will be four output files
1322 for (unsigned i = 0; i < 4; i++)
1323 {
1324 output_stream.push_back(new std::ofstream);
1325 }
1326
1327 // If we are documenting the output, open the files
1328 if (doc_info.is_doc_enabled())
1329 {
1330 std::ostringstream fullname;
1331 fullname << doc_info.directory() << "/hang_nodes_s" << doc_info.number()
1332 << ".dat";
1333 output_stream[0]->open(fullname.str().c_str());
1334 fullname.str("");
1335 fullname << doc_info.directory() << "/hang_nodes_n" << doc_info.number()
1336 << ".dat";
1337 output_stream[1]->open(fullname.str().c_str());
1338 fullname.str("");
1339 fullname << doc_info.directory() << "/hang_nodes_w" << doc_info.number()
1340 << ".dat";
1341 output_stream[2]->open(fullname.str().c_str());
1342 fullname.str("");
1343 fullname << doc_info.directory() << "/hang_nodes_e" << doc_info.number()
1344 << ".dat";
1345 output_stream[3]->open(fullname.str().c_str());
1346 }
1347 }
1348
1349 //================================================================
1350 /// Self test: Check neighbour finding routine. For each element
1351 /// in the tree and for each vertex, determine the
1352 /// distance between the vertex and its position in the
1353 /// neigbour. If the difference is less than
1354 /// Tree::Max_neighbour_finding_tolerance.
1355 /// return success (0), otherwise failure (1)
1356 //=================================================================
1358 {
1359 // Stick pointers to all nodes into Vector and number elements
1360 // in the process
1361 Vector<Tree*> all_forest_nodes_pt;
1362 stick_all_tree_nodes_into_vector(all_forest_nodes_pt);
1363 long int count = 0;
1364 unsigned long num_nodes = all_forest_nodes_pt.size();
1365 for (unsigned long i = 0; i < num_nodes; i++)
1366 {
1367 all_forest_nodes_pt[i]->object_pt()->set_number(++count);
1368 }
1369
1370 // Check neighbours (distance between hanging nodes) -- don't print
1371 // (keep output streams closed)
1372 double max_error = 0.0;
1373 std::ofstream neighbours_file;
1374 std::ofstream neighbours_txt_file;
1376 all_forest_nodes_pt, neighbours_file, neighbours_txt_file, max_error);
1378 {
1379 oomph_info << "\n \n Failed self_test() for QuadTree: Max. error "
1380 << max_error << std::endl
1381 << std::endl;
1382 return 1;
1383 }
1384 else
1385 {
1386 oomph_info << "\n \n Passed self_test() for QuadTree: Max. error "
1387 << max_error << std::endl
1388 << std::endl;
1389 return 0;
1390 }
1391 }
1392
1393 //=================================================================
1394 /// Doc/check all neighbours of quadtree ("nodes") contained in the
1395 /// Vector forest_node_pt. Output into neighbours_file which can
1396 /// be viewed from tecplot with QuadTreeNeighbours.mcr
1397 /// Neighbour info and errors are displayed on
1398 /// neighbours_txt_file. Finally, compute the max. error between
1399 /// vertices when viewed from neighhbouring element.
1400 /// Output is suppressed if the output streams are closed.
1401 //=================================================================
1403 std::ofstream& neighbours_file,
1404 std::ofstream& neighbours_txt_file,
1405 double& max_error)
1406 {
1407 using namespace QuadTreeNames;
1408
1409 int diff_level;
1410 double s_diff;
1411 bool in_neighbouring_tree;
1412 int edge = OMEGA;
1413
1414 Vector<double> s(2);
1415 Vector<double> x(2);
1416 Vector<int> prev_son_type;
1417
1418 Vector<unsigned> translate_s(2);
1419 Vector<double> s_lo(2);
1420 Vector<double> s_hi(2);
1421
1422 Vector<double> x_small(2);
1423 Vector<double> x_large(2);
1424
1425 // Initialise error in vertex positions
1426 max_error = 0.0;
1427
1428
1429 // Loop over all elements to assign numbers for plotting
1430 // -----------------------------------------------------
1431 unsigned long num_nodes = forest_nodes_pt.size();
1432 for (unsigned long i = 0; i < num_nodes; i++)
1433 {
1434 // Set number
1435 forest_nodes_pt[i]->object_pt()->set_number(i);
1436 }
1437
1438 // Loop over all elements for checks
1439 // ---------------------------------
1440 for (unsigned long i = 0; i < num_nodes; i++)
1441 {
1442 // Doc the element itself
1443 QuadTree* el_pt = dynamic_cast<QuadTree*>(forest_nodes_pt[i]);
1444
1445 // If the object is incomplete complain
1446 if (el_pt->object_pt()->nodes_built())
1447 {
1448 // Print it
1449 if (neighbours_file.is_open())
1450 {
1451 dynamic_cast<RefineableQElement<2>*>(el_pt->object_pt())
1452 ->output_corners(neighbours_file, "BLACK");
1453 }
1454
1455 // Loop over directions to find neighbours
1456 // ----------------------------------------
1457 for (int direction = N; direction <= W; direction++)
1458 {
1459 // Initialise difference in levels and coordinate offset
1460 diff_level = 0;
1461 s_diff = 0.0;
1462
1463 // Find greater-or-equal-sized neighbour...
1464 QuadTree* neighb_pt =
1465 el_pt->gteq_edge_neighbour(direction,
1466 translate_s,
1467 s_lo,
1468 s_hi,
1469 edge,
1470 diff_level,
1471 in_neighbouring_tree);
1472
1473 // If neighbour exist and nodes are created: Doc it
1474 if ((neighb_pt != 0) && (neighb_pt->object_pt()->nodes_built()))
1475 {
1476 // Doc neighbour stats
1477 if (neighbours_txt_file.is_open())
1478 {
1479 neighbours_txt_file
1480 << Direct_string[direction] << " neighbour of "
1481 << el_pt->object_pt()->number() << " is "
1482 << neighb_pt->object_pt()->number() << " diff_level "
1483 << diff_level << " s_diff " << s_diff
1484 << " inside neighbour the edge is " << Direct_string[edge]
1485 << std::endl
1486 << std::endl;
1487 }
1488
1489 // Plot neighbour in the appropriate Colour
1490 if (neighbours_file.is_open())
1491 {
1492 dynamic_cast<RefineableQElement<2>*>(neighb_pt->object_pt())
1493 ->output_corners(neighbours_file, Colour[direction]);
1494 }
1495
1496 // Check that local
1497 // coordinates in the larger element (obtained via s_diff)
1498 // lead to the same spatial point as the node vertices
1499 // in the current element
1500 {
1501 if (neighbours_file.is_open())
1502 {
1503 neighbours_file << "ZONE I=2 \n";
1504 }
1505
1506 // (Lower)/(left) vertex:
1507 //-----------------------
1508
1509 // Get coordinates in large (neighbour) element
1510 s[0] = s_lo[0];
1511 s[1] = s_lo[1];
1512 neighb_pt->object_pt()->get_x(s, x_large);
1513
1514 // Get coordinates in small element
1515 Vector<double> s(2);
1516 s[0] = S_base(0, direction);
1517 s[1] = S_base(1, direction);
1518 el_pt->object_pt()->get_x(s, x_small);
1519
1520
1521 // Need to exclude periodic nodes from this check
1522 // There can only be periodic nodes if we have moved into the
1523 // neighbour
1524 bool is_periodic = false;
1525 if (in_neighbouring_tree)
1526 {
1527 // is the node periodic
1528 is_periodic =
1529 el_pt->root_pt()->is_neighbour_periodic(direction);
1530 }
1531
1532 double error = 0.0;
1533 // Only bother to calculate the error if the node is NOT periodic
1534 if (is_periodic == false)
1535 {
1536 for (int i = 0; i < 2; i++)
1537 {
1538 error += pow(x_small[i] - x_large[i], 2);
1539 }
1540 }
1541
1542 // Take the root of the square error
1543 error = sqrt(error);
1544 if (neighbours_txt_file.is_open())
1545 {
1546 neighbours_txt_file << "Error (1) " << error << std::endl;
1547 }
1548
1549 if (std::fabs(error) > max_error)
1550 {
1551 max_error = std::fabs(error);
1552 }
1553
1554 if (neighbours_file.is_open())
1555 {
1556 neighbours_file << x_large[0] << " " << x_large[1] << " 0 \n";
1557 }
1558
1559 // (Upper)/(right) vertex:
1560 //------------------------
1561
1562 // Get coordinates in large (neighbour) element
1563 s[0] = s_hi[0];
1564 s[1] = s_hi[1];
1565 neighb_pt->object_pt()->get_x(s, x_large);
1566
1567 // Get coordinates in small element
1568 s[0] = S_base(0, direction) + S_step(0, direction);
1569 s[1] = S_base(1, direction) + S_step(1, direction);
1570 el_pt->object_pt()->get_x(s, x_small);
1571
1572 error = 0.0;
1573 // Only do this if we are NOT periodic
1574 if (is_periodic == false)
1575 {
1576 for (int i = 0; i < 2; i++)
1577 {
1578 error += pow(x_small[i] - x_large[i], 2);
1579 }
1580 }
1581 // Take the root of the square error
1582 error = sqrt(error);
1583
1584 // error=
1585 // sqrt(pow(x_small[0]-x_large[0],2)+pow(x_small[1]-x_large[1],2));
1586 if (neighbours_txt_file.is_open())
1587 {
1588 neighbours_txt_file << "Error (2) " << error << std::endl;
1589 }
1590
1591 if (std::fabs(error) > max_error)
1592 {
1593 max_error = std::fabs(error);
1594 }
1595
1596 if (neighbours_file.is_open())
1597 {
1598 neighbours_file << x_large[0] << " " << x_large[1] << " 0 \n";
1599 }
1600 }
1601 // else
1602 // {
1603 // // No neighbour: write dummy zone so tecplot can find
1604 // four
1605 // // neighbours for every element
1606 // if (neighbours_file.is_open())
1607 // {
1608 // neighbours_file << "ZONE I=2 \n";
1609 // neighbours_file << "-0.05 -0.05 0 \n";
1610 // neighbours_file << "-0.05 -0.05 0 \n";
1611 // }
1612 // }
1613 }
1614 // If neighbour does not exist: Insert blank zones into file
1615 // so that tecplot can find four neighbours for every element
1616 else
1617 {
1618 if (neighbours_file.is_open())
1619 {
1620 neighbours_file << "ZONE \n 0.00 0.00 0 \n";
1621 neighbours_file << "ZONE I=2 \n";
1622 neighbours_file << "-0.05 -0.05 0 \n";
1623 neighbours_file << "-0.05 -0.05 0 \n";
1624 }
1625 }
1626 }
1627 } // End of case when element can be documented
1628 }
1629 }
1630
1631} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
Information for documentation of results: Directory and file number to enable output in the form RESL...
bool is_doc_enabled() const
Are we documenting?
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1885
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
Base class for all quad elements.
Definition: Qelements.h:821
QuadTreeForest()
Default constructor (empty and broken)
Definition: quadtree.h:412
QuadTreeRoot * quadtree_pt(const unsigned &i)
Return pointer to i-th root quadtree in this forest. (Performs a dynamic cast from the TreeRoot to a ...
Definition: quadtree.h:462
void check_all_neighbours(DocInfo &doc_info)
Document and check all the neighbours of all the nodes in the forest. DocInfo object specifies the ou...
Definition: quadtree.cc:1245
void find_neighbours()
Construct the neighbour lookup scheme.
Definition: quadtree.cc:891
unsigned self_test()
Self-test: Check all neighbours. Return success (0) if the max. distance between corresponding points...
Definition: quadtree.cc:1357
void construct_north_equivalents()
Construct the rotation schemes.
Definition: quadtree.cc:1067
void open_hanging_node_files(DocInfo &doc_info, Vector< std::ofstream * > &output_stream)
Open output files that will store any hanging nodes in the forest and return a vector of the streams.
Definition: quadtree.cc:1318
QuadTreeRoot * quad_neigh_pt(const unsigned &i, const int &direction)
Given the number i of the root quadtree in this forest, return pointer to its neighbour in the specif...
Definition: quadtree.h:471
QuadTreeRoot is a QuadTree that forms the root of a (recursive) quadtree. The "root node" is special ...
Definition: quadtree.h:293
int & north_equivalent(const int &neighbour)
Return north equivalent of the neighbours in specified direction: When viewed from the current quadtr...
Definition: quadtree.h:354
int direction_of_neighbour(QuadTreeRoot *quadtree_root_pt)
If quadtree_root_pt is a neighbour, return the direction [N/S/E/W] in which it is found,...
Definition: quadtree.h:377
QuadTree class: Recursively defined, generalised quadtree.
Definition: quadtree.h:104
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: quadtree.h:199
static void doc_neighbours(Vector< Tree * > forest_nodes_pt, std::ofstream &neighbours_file, std::ofstream &neighbours_txt_file, double &max_error)
Doc/check all neighbours of quadtree (nodes) contained in the Vector forest_node_pt....
Definition: quadtree.cc:1402
static DenseMatrix< double > S_step
S_step(i,direction) Increments for coordinate s[i] when progressing along the edge indicated by direc...
Definition: quadtree.h:255
static Vector< int > Reflect_edge
Get opposite edge, e.g. Reflect_edge[N]=S.
Definition: quadtree.h:258
static DenseMatrix< int > Reflect
Reflection scheme: Reflect(direction,quadrant): Get mirror of quadrant in specified direction....
Definition: quadtree.h:267
static DenseMatrix< int > Rotate
Rotate coordinates: If North becomes NorthIs then direction becomes Rotate(NorthIs,...
Definition: quadtree.h:271
static void setup_static_data()
Setup the static data, rotation and reflection schemes, etc.
Definition: quadtree.cc:120
static DenseMatrix< double > S_base
S_base(i,direction): Initial value for coordinate s[i] on the edge indicated by direction (S/E/N/W)
Definition: quadtree.h:250
static DenseMatrix< int > S_direct
S_direct(direction,son_quadrant): The lower left corner of son_quadrant has an offset of h/2 S_direct...
Definition: quadtree.h:282
unsigned self_test()
Self-test: Check all neighbours. Return success (0) if the max. distance between corresponding points...
Definition: quadtree.cc:814
static Vector< std::string > Colour
Colours for neighbours in various directions.
Definition: quadtree.h:246
QuadTree * gteq_edge_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level, bool &in_neighbouring_tree) const
Return pointer to greater or equal-sized edge neighbour in specified direction; also provide info reg...
Definition: quadtree.cc:413
static DenseMatrix< bool > Is_adjacent
Array of direction/quadrant adjacency scheme: Is_adjacent(i_vertex_or_edge,j_quadrant): Is edge/verte...
Definition: quadtree.h:263
void stick_neighbouring_leaves_into_vector(Vector< const QuadTree * > &tree_neighbouring_nodes, Vector< Vector< double > > &tree_neighbouring_s_lo, Vector< Vector< double > > &tree_neighbouring_s_hi, Vector< int > &tree_neighbouring_diff_level, const QuadTree *my_neigh_pt, const int &direction) const
Traverse Tree: Preorder traverse and stick pointers to neighbouring leaf nodes (only) into Vector.
Definition: quadtree.cc:751
static bool Static_data_has_been_setup
Bool indicating that static member data has been setup.
Definition: quadtree.h:232
static DenseMatrix< int > Rotate_angle
Angle betwen rotated coordinates: If old_direction becomes new_direction then the angle between the a...
Definition: quadtree.h:277
long number() const
Element number (for debugging/plotting)
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
Refineable version of QElement<2,NNODE_1D>.
A TreeForest consists of a collection of TreeRoots. Each member tree can have neighbours in various e...
Definition: tree.h:409
void stick_all_tree_nodes_into_vector(Vector< Tree * > &all_forest_nodes)
Traverse forest and stick pointers to all "nodes" into Vector.
Definition: tree.cc:405
unsigned ntree()
Number of trees in forest.
Definition: tree.h:460
Vector< TreeRoot * > Trees_pt
Vector containing the pointers to the trees.
Definition: tree.h:480
bool is_neighbour_periodic(const int &direction)
Return whether the neighbour in the particular direction is periodic.
Definition: tree.h:364
TreeRoot *& neighbour_pt(const int &direction)
Return the pointer to the neighbouring TreeRoots in specified direction. Returns NULL if there's no n...
Definition: tree.h:357
Tree * Father_pt
Pointer to the Father of the Tree.
Definition: tree.h:296
void stick_all_tree_nodes_into_vector(Vector< Tree * > &)
Traverse and stick pointers to all "nodes" into Vector.
Definition: tree.cc:277
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:141
TreeRoot * Root_pt
Pointer to the root of the tree.
Definition: tree.h:292
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:88
int Son_type
Son type (e.g. SW/SE/NW/NE in a quadtree)
Definition: tree.h:305
static double & max_neighbour_finding_tolerance()
Max. allowed discrepancy in neighbour finding routine (distance between points when identified from t...
Definition: tree.h:255
int Level
Level of the Tree (level 0 = root)
Definition: tree.h:302
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
Vector< Tree * > Son_pt
Vector of pointers to the sons of the Tree.
Definition: tree.h:299
static double Max_neighbour_finding_tolerance
Max. allowed discrepancy in neighbour finding routine (distance between points when identified from t...
Definition: tree.h:313
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...