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