binary_tree.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 BinaryTree and BinaryTreeForest
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 "binary_tree.h"
43 
44 namespace oomph
45 {
46  //========================================================================
47  /// Boolean indicating that static member data has been setup.
48  //========================================================================
50 
51  //========================================================================
52  /// Colours for neighbours in various directions (static data).
53  //========================================================================
54  Vector<std::string> BinaryTree::Colour;
55 
56  //========================================================================
57  /// S_base(direction): Initial value for coordinate s on the edge
58  /// indicated by direction (L/R) (static data).
59  //========================================================================
60  Vector<double> BinaryTree::S_base;
61 
62  //========================================================================
63  /// Translate (enumerated) directions into strings (static data).
64  //========================================================================
65  Vector<std::string> BinaryTree::Direct_string;
66 
67  //========================================================================
68  /// Get opposite edge, e.g. Reflect_edge[N]=S (static data)
69  //========================================================================
70  Vector<int> BinaryTree::Reflect_edge;
71 
72  //========================================================================
73  /// Array of direction/segment adjacency scheme:
74  /// Is_adjacent(i_vertex,j_segment): Is vertex adjacent to segment?
75  /// (static data)
76  //========================================================================
77  DenseMatrix<bool> BinaryTree::Is_adjacent;
78 
79  //========================================================================
80  /// Reflection scheme: Reflect(direction,segment): Get mirror of segment
81  /// in specified direction. E.g. Reflect(L,L)=R (static data)
82  //========================================================================
83  DenseMatrix<int> BinaryTree::Reflect;
84 
85  //========================================================================
86  /// Set up the static data stored in the BinaryTree -- this needs to be
87  /// called before BinaryTrees can be used. Automatically called by
88  /// RefineableLineMesh constructor.
89  //========================================================================
91  {
92  using namespace BinaryTreeNames;
93 
94 #ifdef PARANOID
96  {
97  std::ostringstream error_stream;
98  error_stream
99  << "Inconsistent enumeration! \n Tree::OMEGA=" << Tree::OMEGA
100  << "\nBinaryTree::OMEGA=" << BinaryTree::OMEGA
101  << std::endl;
102  throw OomphLibError(error_stream.str(),
103  OOMPH_CURRENT_FUNCTION,
104  OOMPH_EXCEPTION_LOCATION);
105  }
106 #endif
107 
108  // Set flag to indicate that static data has been setup
110 
111  // Tecplot colours for neighbours in various directions
112  Colour.resize(27);
113 
114  Colour[L]="RED";
115  Colour[R]="CYAN";
116  Colour[OMEGA]="YELLOW";
117 
118  // S_base(direction): Initial value for coordinate s on the
119  // edge indicated by direction (L/R)
120  S_base.resize(27);
121 
122  S_base[L]=-1.0;
123  S_base[R]=1.0;
124 
125  // Translate (enumerated) directions into strings
126  Direct_string.resize(27);
127 
128  Direct_string[L]="L";
129  Direct_string[R]="R";
130  Direct_string[OMEGA]="OMEGA";
131 
132  // Build direction/segment adjacency scheme: Is_adjacent(i_vertex,j_segment):
133  // Is vertex adjacent to segment?
134  Is_adjacent.resize(27,27);
135 
136  Is_adjacent(L,L)=true;
137  Is_adjacent(R,L)=false;
138  Is_adjacent(L,R)=false;
139  Is_adjacent(R,R)=true;
140 
141  // Build reflection scheme: Reflect(direction,segment):
142  // Get mirror of segment in direction
143  Reflect.resize(27,27);
144 
145  Reflect(L,L)=R;
146  Reflect(R,L)=R;
147  Reflect(L,R)=L;
148  Reflect(R,R)=L;
149 
150  // Get opposite edge, e.g. Reflect_edge(L)=R
151  Reflect_edge.resize(27);
152 
153  Reflect_edge[L]=R;
154  Reflect_edge[R]=L;
155  }
156 
157  //========================================================================
158  /// Return pointer to greater or equal-sized edge neighbour in specified
159  /// \c direction; also provide info regarding the relative size of the
160  /// neighbour:
161  /// - In the present binary tree, the left vertex is located at the local
162  /// coordinate s = -1. This point is located at the local coordinate
163  /// s = \c s_in_neighbour[0] in the neighbouring binary tree.
164  /// - We're looking for a neighbour in the specified \c direction. When
165  /// viewed from the neighbouring binary tree, the edge that separates
166  /// the present binary tree from its neighbour is the neighbour's
167  /// \c edge edge. Since in 1D there can be no rotation between the two
168  /// binary trees, this is a simple reflection. For instance, if we're
169  /// looking for a neighhbour in the \c L [eft] \c direction, \c edge
170  /// will be \c R [ight].
171  /// - \c diff_level <= 0 indicates the difference in refinement levels
172  /// between the two neighbours. If \c diff_level==0, the neighbour has
173  /// the same size as the current binary tree.
174  /// - \c in_neighbouring_tree returns true if we have had to flip to a
175  /// different root, even if that root is actually the same (as it can
176  /// be in periodic problems).
177  //========================================================================
179  gteq_edge_neighbour(const int& direction,
180  Vector<double>& s_in_neighbour,
181  int& edge, int& diff_level,
182  bool &in_neighbouring_tree) const
183  {
184  using namespace BinaryTreeNames;
185 
186 #ifdef PARANOID
187  if ((direction!=L)&&(direction!=R))
188  {
189  std::ostringstream error_stream;
190  error_stream << "Direction " << direction
191  << " is not L or R" << std::endl;
192 
193  throw OomphLibError(error_stream.str(),
194  OOMPH_CURRENT_FUNCTION,
195  OOMPH_EXCEPTION_LOCATION);
196  }
197 #endif
198 
199  // Initialise in_neighbouring tree to false. It will be set to true
200  // during the recursion if we do actually hop over into the neighbour
201  in_neighbouring_tree = false;
202 
203  // Maximum level to which we're allowed to descend (we only want
204  // greater-or-equal-sized neighbours)
205  int max_level = Level;
206 
207  // Current element has the following root:
208  BinaryTreeRoot* orig_root_pt = dynamic_cast<BinaryTreeRoot*>(Root_pt);
209 
210  // Initialise offset in local coordinate
211  double s_diff = 0;
212 
213  // Initialise difference in level
214  diff_level = 0;
215 
216  // Find neighbour
217  BinaryTree* return_pt = gteq_edge_neighbour(direction,s_diff,diff_level,
218  in_neighbouring_tree,
219  max_level,orig_root_pt);
220 
221  BinaryTree* neighb_pt = return_pt;
222 
223  // If neighbour exists...
224  if (neighb_pt!=0)
225  {
226  // What's the direction of the interfacial edge when viewed from within
227  // the neighbour element?
228  edge = Reflect_edge[direction];
229 
230  // What's the local coordinate s at which this edge is located in the
231  // neighbouring element?
232  s_in_neighbour[0] = S_base[edge];
233  }
234  return return_pt;
235  }
236 
237  //========================================================================
238  /// Find `greater-or-equal-sized edge neighbour' in given direction (L/R).
239  ///
240  /// This is an auxiliary routine which allows neighbour finding in
241  /// adjacent binary trees. Needs to keep track of previous son types
242  /// and the maximum level to which search is performed.
243  ///
244  /// Parameters:
245  /// - direction (L/R): Direction in which neighbour has to be found.
246  /// - s_diff: Offset of left vertex from corresponding vertex in
247  /// neighbour. Note that this is input/output as it needs to be
248  /// incremented/decremented during the recursive calls to this function.
249  /// - edge: We're looking for the neighbour across our edge 'direction'
250  /// (L/R). When viewed from the neighbour, this edge is `edge' (L/R).
251  /// Since there is no relative rotation between neighbours this is a
252  /// mere reflection, e.g. direction=L --> edge=R etc.
253  /// - diff_level <= 0 indicates the difference in binary tree levels
254  /// between the current element and its neighbour.
255  /// - max_level is the maximum level to which the neighbour search is
256  /// allowed to proceed. This is again necessary because in a forest,
257  /// the neighbour search isn't based on pure recursion.
258  /// - orig_root_pt identifies the root node of the element whose
259  /// neighbour we're really trying to find by all these recursive calls.
260  //========================================================================
262  const int& direction,
263  double& s_diff,
264  int& diff_level,
265  bool &in_neighbouring_tree,
266  int max_level,
267  BinaryTreeRoot* const &orig_root_pt) const
268  {
269  using namespace BinaryTreeNames;
270 
271 #ifdef PARANOID
272  if ((direction!=L)&&(direction!=R))
273  {
274  std::ostringstream error_stream;
275  error_stream << "Direction " << direction
276  << " is not L or R" << std::endl;
277 
278  throw OomphLibError(error_stream.str(),
279  OOMPH_CURRENT_FUNCTION,
280  OOMPH_EXCEPTION_LOCATION);
281  }
282 #endif
283 
284  BinaryTree* next_el_pt;
285  BinaryTree* return_el_pt;
286 
287  // STEP 1: Find the neighbour's father
288  // -------
289 
290  // Does the element have a father?
291  if(Father_pt!=0)
292  {
293  // If the present segment (whose location inside its father element
294  // is specified by Son_type) is adjacent to the father's edge in the
295  // required direction, then its neighbour has a different father
296  // ---> we need to climb up the tree to the father and find his
297  // neighbour in the required direction. Note that this is the cunning
298  // recursive part. The returning may not stop until we hit the very
299  // top of the tree, when the element does NOT have a father.
300  if(Is_adjacent(direction,Son_type))
301  {
302  next_el_pt = dynamic_cast<BinaryTree*>(Father_pt)
303  ->gteq_edge_neighbour(direction,s_diff,diff_level,
304  in_neighbouring_tree,
305  max_level,orig_root_pt);
306  }
307  // If the present segment is not adjacent to the father's edge in the
308  // required direction, then the neighbour has the same father and is
309  // obtained by the appropriate reflection inside the father element.
310  // This will only be called if we have not left the original tree.
311  else
312  {
313  next_el_pt = dynamic_cast<BinaryTree*>(Father_pt);
314  }
315 
316  // We're about to ascend one level
317  diff_level -= 1;
318 
319  // Work out position of left corner of present edge in its father element
320  s_diff += pow(0.5,-diff_level);
321 
322  // STEP 2: We have now located the neighbour's father and need to
323  // ------- find the appropriate son.
324 
325  // Buffer cases where the neighbour (and hence its father) lie outside
326  // the boundary
327  if(next_el_pt!=0)
328  {
329  // If the father is a leaf then we can't descend to the same
330  // level as the present node ---> simply return the father himself
331  // as the (greater) neighbour. Same applies if we are about to
332  // descend lower than the max_level (in a neighbouring tree).
333  if((next_el_pt->Son_pt.size()==0)||(next_el_pt->Level>max_level-1))
334  {
335  return_el_pt = next_el_pt;
336  }
337  // We have located the neighbour's father: The position of the
338  // neighbour is obtained by `reflecting' the position of the
339  // node itself. We know exactly how to reflect because we know which
340  // son type we are and we have the pointer to the neighbour's father.
341  else
342  {
343  int son_segment = Reflect(direction,Son_type);
344 
345  // The next element in the tree is the appropriate son of the
346  // neighbour's father
347  return_el_pt
348  = dynamic_cast<BinaryTree*>(next_el_pt->Son_pt[son_segment]);
349 
350  // Work out position of left corner of present edge in next
351  // higher element
352  s_diff -= pow(0.5,-diff_level);
353 
354  // We have just descended one level
355  diff_level += 1;
356  }
357  }
358  // The neighbour's father lies outside the boundary --> the neighbour
359  // itself does too --> return NULL
360  else
361  {
362  return_el_pt=0;
363  }
364  }
365  // Element does not have a father --> check if it has a neighbouring
366  // tree in the appropriate direction
367  else
368  {
369  // Find neighbouring root
370  if(Root_pt->neighbour_pt(direction)!=0)
371  {
372  // In this case we have moved to a neighbour, so set the flag
373  in_neighbouring_tree = true;
374  return_el_pt
375  = dynamic_cast<BinaryTreeRoot*>(Root_pt->neighbour_pt(direction));
376  }
377  // No neighbouring tree, so there really is no neighbour --> return NULL
378  else
379  {
380  return_el_pt=0;
381  }
382  }
383 
384  return return_el_pt;
385  }
386 
387  //========================================================================
388  /// Self-test: Check neighbour finding routine. For each element in the
389  /// tree and for each vertex, determine the distance between the vertex
390  /// and its position in the neighbour. If the difference is less than
391  /// Tree::Max_neighbour_finding_tolerance return success (0), otherwise
392  /// failure (1).
393  //========================================================================
395  {
396  // Stick pointers to all nodes into Vector and number elements
397  // in the process
398  Vector<Tree*> all_nodes_pt;
399  stick_all_tree_nodes_into_vector(all_nodes_pt);
400  long int count=0;
401  unsigned long num_nodes=all_nodes_pt.size();
402  for (unsigned long i=0;i<num_nodes;i++)
403  {
404  all_nodes_pt[i]->object_pt()->set_number(++count);
405  }
406 
407  // Check neighbours (distance between hanging nodes) -- don't print
408  // (keep output streams closed)
409  double max_error=0.0;
410  std::ofstream neighbours_file;
411  std::ofstream neighbours_txt_file;
412  BinaryTree::doc_neighbours(all_nodes_pt,neighbours_file,
413  neighbours_txt_file, max_error);
414 
415  if (max_error>Max_neighbour_finding_tolerance)
416  {
417  oomph_info << "\n \n Failed self_test() for BinaryTree: Max. error "
418  << max_error << std::endl<< std::endl;
419  return 1;
420  }
421  else
422  {
423  oomph_info << "\n \n Passed self_test() for BinaryTree: Max. error "
424  << max_error << std::endl<< std::endl;
425  return 0;
426  }
427  }
428 
429 
430 
431 
432 //========================================================================
433 /// Constructor for BinaryTreeForest:
434 ///
435 /// Pass:
436 /// - trees_pt[], the Vector of pointers to the constituent trees
437 /// (BinaryTreeRoot objects)
438 //========================================================================
440  TreeForest(trees_pt)
441  {
442 #ifdef LEAK_CHECK
443  LeakCheckNames::BinaryTreeForest_build+=1;
444 #endif
445 
446  using namespace BinaryTreeNames;
447 
448  // Set up the neighbours
449  find_neighbours();
450  }
451 
452  //========================================================================
453  /// Set up the neighbour lookup schemes for all constituent binary trees.
454  //========================================================================
456  {
457  using namespace BinaryTreeNames;
458 
459  unsigned numtrees = ntree();
460  unsigned n = 0; // to store nnode1d
461  if(numtrees>0)
462  {
463  n=Trees_pt[0]->object_pt()->nnode_1d();
464  }
465  else
466  {
467  throw OomphLibError(
468  "Trying to setup the neighbour scheme for an empty forest\n",
469  OOMPH_CURRENT_FUNCTION,
470  OOMPH_EXCEPTION_LOCATION);
471  }
472 
473  // Number of vertex nodes: 2
474  unsigned n_vertex_node = 2;
475 
476  // Find connected trees by identifying those whose associated elements
477  // share a common vertex node
478  std::map<Node*,std::set<unsigned> > tree_assoc_with_vertex_node;
479 
480  // Loop over all trees
481  for(unsigned i=0;i<numtrees;i++)
482  {
483  // Loop over the vertex nodes of the associated element
484  for (unsigned j=0;j<n_vertex_node;j++)
485  {
486  Node* nod_pt = dynamic_cast< LineElementBase*>
487  (Trees_pt[i]->object_pt())->vertex_node_pt(j);
488  tree_assoc_with_vertex_node[nod_pt].insert(i);
489  }
490  }
491 
492  // For each tree we store a set of neighbouring trees
493  // i.e. trees that share a node
494  Vector<std::set<unsigned> > neighb_tree(numtrees);
495 
496  // Loop over vertex nodes
497  for (std::map<Node*,std::set<unsigned> >::iterator it=
498  tree_assoc_with_vertex_node.begin();
499  it!=tree_assoc_with_vertex_node.end();it++)
500  {
501  // Loop over connected elements twice
502  for (std::set<unsigned>::iterator it_el1=it->second.begin();
503  it_el1!=it->second.end();it_el1++)
504  {
505  unsigned i=(*it_el1);
506  for (std::set<unsigned>::iterator it_el2=it->second.begin();
507  it_el2!=it->second.end();it_el2++)
508  {
509  unsigned j=(*it_el2);
510  // These two elements are connected
511  if (i!=j)
512  {
513  neighb_tree[i].insert(j);
514  }
515  }
516  }
517  }
518 
519  // Loop over all trees
520  for(unsigned i=0;i<numtrees;i++)
521  {
522  // Loop over their neighbours
523  for(std::set<unsigned>::iterator it=neighb_tree[i].begin();
524  it!=neighb_tree[i].end();it++)
525  {
526  unsigned j=(*it);
527 
528  // is it the left-hand neighbour?
529  bool is_L_neighbour=
530  Trees_pt[j]->object_pt()->get_node_number(
531  Trees_pt[i]->object_pt()->node_pt(0))!=-1;
532 
533  // is it the right-hand neighbour?
534  bool is_R_neighbour=
535  Trees_pt[j]->object_pt()->get_node_number(
536  Trees_pt[i]->object_pt()->node_pt(n-1))!=-1;
537 
538  if(is_L_neighbour) Trees_pt[i]->neighbour_pt(L)=Trees_pt[j];
539  if(is_R_neighbour) Trees_pt[i]->neighbour_pt(R)=Trees_pt[j];
540  }
541  } // End of loop over all trees
542  }
543 
544  //========================================================================
545  /// Document and check all the neighbours in all the nodes in the forest.
546  //========================================================================
548  {
549  // Create Vector of elements
550  Vector<Tree*> all_tree_nodes_pt;
551  this->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
552 
553  // Create storage for information files
554  std::ofstream neigh_file;
555  std::ofstream neigh_txt_file;
556 
557  // If we are documenting the results, then open the files
558  if(doc_info.is_doc_enabled())
559  {
560  std::ostringstream fullname;
561  fullname << doc_info.directory() << "/neighbours"
562  << doc_info.number() << ".dat";
563  oomph_info << "opened " << fullname.str() << " to doc neighbours"
564  << std::endl;
565  neigh_file.open(fullname.str().c_str());
566  fullname.str("");
567  fullname << doc_info.directory() << "/neighbours"
568  << doc_info.number() << ".txt";
569  oomph_info << "opened " << fullname.str() << " to doc neighbours"
570  << std::endl;
571  neigh_txt_file.open(fullname.str().c_str());
572  }
573 
574  // Call the standard documentation function
575  double max_error=0.0;
576  BinaryTree::doc_neighbours(all_tree_nodes_pt,neigh_file,
577  neigh_txt_file,max_error);
578 
579  // If the error is too large, complain
581  {
582  std::ostringstream error_stream;
583  error_stream << "Max. error in binary tree neighbour finding: "
584  << max_error << " is too big" << std::endl;
585  error_stream
586  << "i.e. bigger than Tree::max_neighbour_finding_tolerance()="
587  << Tree::max_neighbour_finding_tolerance() << std::endl;
588 
589  // Close the files if they were opened
590  if(doc_info.is_doc_enabled())
591  {
592  neigh_file.close();
593  neigh_txt_file.close();
594  }
595 
596  throw OomphLibError(error_stream.str(),
597  OOMPH_CURRENT_FUNCTION,
598  OOMPH_EXCEPTION_LOCATION);
599  }
600  else
601  {
602  oomph_info << "Max. error in binary tree neighbour finding: "
603  << max_error << " is OK" << std::endl;
604  oomph_info
605  << "i.e. less than BinaryTree::max_neighbour_finding_tolerance()="
607  }
608 
609  // Close the files if they were opened
610  if(doc_info.is_doc_enabled())
611  {
612  neigh_file.close();
613  neigh_txt_file.close();
614  }
615  }
616 
617  //========================================================================
618  /// Self test: Check neighbour finding routine. For each element in the
619  /// tree and for each vertex, determine the distance between the vertex
620  /// and its position in the neighbour. If the difference is less than
621  /// Tree::Max_neighbour_finding_tolerance return success (0), otherwise
622  /// failure (1).
623  //========================================================================
625  {
626  // Stick pointers to all nodes into Vector and number elements
627  // in the process
628  Vector<Tree*> all_forest_nodes_pt;
629  stick_all_tree_nodes_into_vector(all_forest_nodes_pt);
630  long int count=0;
631  unsigned long num_nodes=all_forest_nodes_pt.size();
632  for(unsigned long i=0;i<num_nodes;i++)
633  {
634  all_forest_nodes_pt[i]->object_pt()->set_number(++count);
635  }
636 
637  // Check neighbours (distance between hanging nodes) -- don't print
638  // (keep output streams closed)
639  double max_error=0.0;
640  std::ofstream neighbours_file;
641  std::ofstream neighbours_txt_file;
642  BinaryTree::doc_neighbours(all_forest_nodes_pt,neighbours_file,
643  neighbours_txt_file, max_error);
645  {
646  oomph_info << "\n \n Failed self_test() for BinaryTree: Max. error "
647  << max_error << std::endl<< std::endl;
648  return 1;
649  }
650  else
651  {
652  oomph_info << "\n \n Passed self_test() for BinaryTree: Max. error "
653  << max_error << std::endl<< std::endl;
654  return 0;
655  }
656  }
657 
658  //========================================================================
659  /// Doc/check all neighbours of binary tree ("nodes") contained in the
660  /// Vector forest_node_pt. Output into neighbours_file which can be viewed
661  /// from tecplot with BinaryTreeNeighbours.mcr. Neighbour info and errors
662  /// are displayed on neighbours_txt_file. Finally, compute the maximum
663  /// error between vertices when viewed from neighbouring element. Output
664  /// is suppressed if the output streams are closed.
665  //========================================================================
667  std::ofstream& neighbours_file,
668  std::ofstream& neighbours_txt_file,
669  double& max_error)
670  {
671  using namespace BinaryTreeNames;
672 
673  int diff_level;
674  double s_diff;
675  bool in_neighbouring_tree;
676  int edge=OMEGA;
677 
678  Vector<double> s(1);
679  Vector<double> x(1);
680  Vector<int> prev_son_type;
681 
682  Vector<double> s_in_neighbour(1);
683 
684  Vector<double> x_small(1);
685  Vector<double> x_large(1);
686 
687  // Initialise error in vertex positions
688  max_error=0.0;
689 
690  // Loop over all elements to assign numbers for plotting
691  // -----------------------------------------------------
692  unsigned long num_nodes=forest_nodes_pt.size();
693  for(unsigned long i=0;i<num_nodes;i++)
694  {
695  // Set number
696  forest_nodes_pt[i]->object_pt()->set_number(i);
697  }
698 
699  // Loop over all elements for checks
700  // ---------------------------------
701  for(unsigned long i=0;i<num_nodes;i++)
702  {
703 
704  // Doc the element itself
705  BinaryTree* el_pt=dynamic_cast<BinaryTree*>(forest_nodes_pt[i]);
706 
707  // If the object is incomplete, complain
708  if(el_pt->object_pt()->nodes_built())
709  {
710  // Print it
711  if (neighbours_file.is_open())
712  {
713  dynamic_cast<RefineableQElement<1>*>(el_pt
714  ->object_pt())->output_corners(neighbours_file,"BLACK");
715  }
716 
717  // Loop over directions to find neighbours
718  // ---------------------------------------
719  for (int direction=L;direction<=R;direction++)
720  {
721  // Initialise difference in levels and coordinate offset
722  diff_level=0;
723  s_diff=0.0;
724 
725  // Find greater-or-equal-sized neighbour...
726  BinaryTree* neighb_pt=el_pt->gteq_edge_neighbour(direction,
727  s_in_neighbour,
728  edge,diff_level,
729  in_neighbouring_tree);
730 
731  // If neighbour exist and nodes are created: Doc it
732  if((neighb_pt!=0) && (neighb_pt->object_pt()->nodes_built()))
733  {
734  // Doc neighbour stats
735  if(neighbours_txt_file.is_open())
736  {
737  neighbours_txt_file << Direct_string[direction]
738  << " neighbour of " <<
739  el_pt->object_pt()->number() << " is " <<
740  neighb_pt->object_pt()->number()
741  << " diff_level " << diff_level
742  << " s_diff " << s_diff
743  << " inside neighbour the edge is "
744  << Direct_string[edge] << std::endl
745  << std::endl;
746  }
747 
748  // Plot neighbour in the appropriate colour
749  if(neighbours_file.is_open())
750  {
751  dynamic_cast<RefineableQElement<1>*>(neighb_pt->object_pt())->
752  output_corners(neighbours_file,Colour[direction]);
753  }
754 
755  // Check that local coordinates in the larger element (obtained
756  // via s_diff) lead to the same spatial point as the node vertices
757  // in the current element
758  {
759  if(neighbours_file.is_open())
760  {
761  neighbours_file << "ZONE I=1 \n";
762  }
763 
764  // Left vertex:
765  // ------------
766 
767  // Get coordinates in large (neighbour) element
768  s[0] = s_in_neighbour[0];
769  neighb_pt->object_pt()->get_x(s,x_large);
770 
771  // Get coordinates in small element
772  Vector<double> s(1);
773  s[0] = S_base[direction];
774  el_pt->object_pt()->get_x(s,x_small);
775 
776  // Need to exclude periodic nodes from this check. There can
777  // only be periodic nodes if we have moved into the neighbour
778  bool is_periodic = false;
779  if(in_neighbouring_tree)
780  {
781  // Is the node periodic?
782  is_periodic = el_pt->root_pt()->is_neighbour_periodic(direction);
783  }
784 
785  double error = 0.0;
786  // Only bother to calculate the error if the node is NOT periodic
787  if(is_periodic==false)
788  {
789  error += pow(x_small[0]-x_large[0],2);
790  }
791 
792  // Take the root of the square error
793  error = sqrt(error);
794  if(neighbours_txt_file.is_open())
795  {
796  neighbours_txt_file << "Error (1) " << error << std::endl;
797  }
798 
799  if(std::fabs(error)>max_error)
800  {
801  max_error=std::fabs(error);
802  }
803 
804  if(neighbours_file.is_open())
805  {
806  neighbours_file << x_large[0] << " 0 \n";
807  }
808 
809  // Right vertex:
810  // -------------
811 
812  // Get coordinates in large (neighbour) element
813  s[0] = s_in_neighbour[0];
814  neighb_pt->object_pt()->get_x(s,x_large);
815 
816  // Get coordinates in small element
817  s[0] = S_base[direction];
818  el_pt->object_pt()->get_x(s,x_small);
819 
820  error = 0.0;
821  // Only do this if we are NOT periodic
822  if(is_periodic==false)
823  {
824  error += pow(x_small[0]-x_large[0],2);
825  }
826  // Take the root of the square error
827  error = sqrt(error);
828 
829  // error =
830  // sqrt(pow(x_small[0]-x_large[0],2)+pow(x_small[1]-x_large[1],2));
831  if(neighbours_txt_file.is_open())
832  {
833  neighbours_txt_file << "Error (2) " << error << std::endl;
834  }
835 
836  if(std::fabs(error)>max_error)
837  {
838  max_error=std::fabs(error);
839  }
840 
841  if(neighbours_file.is_open())
842  {
843  neighbours_file << x_large[0] << " 0 \n";
844  }
845 
846  }
847 // else
848 // {
849 // // No neighbour: write dummy zone so tecplot can find four
850 // // neighbours for every element
851 // if (neighbours_file.is_open())
852 // {
853 // neighbours_file << "ZONE I=1 \n";
854 // neighbours_file << "-0.05 -0.05 0 \n";
855 // neighbours_file << "-0.05 -0.05 0 \n";
856 // }
857 // }
858  }
859  // If neighbour does not exist: Insert blank zones into file
860  // so that tecplot can find four neighbours for every element
861  else
862  {
863  if (neighbours_file.is_open())
864  {
865  neighbours_file << "ZONE \n 0.00 0 \n";
866  neighbours_file << "ZONE I=1 \n";
867  neighbours_file << "-0.05 0 \n";
868  neighbours_file << "-0.05 0 \n";
869  }
870  }
871  }
872  } // End of case when element can be documented
873  }
874  }
875 
876 } // End of namespace
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
Base class for all line elements.
Definition: Qelements.h:504
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
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
static bool Static_data_has_been_setup
Boolean indicating that static member data has been setup.
Definition: binary_tree.h:193
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 binary tree (nodes) contained in the Vector forest_node_pt. Output into neighbours_file which can be viewed from tecplot with BinaryTreeNeighbours.mcr. Neighbour info and errors are displayed on neighbours_txt_file. Finally, compute the maximum error between vertices when viewed from the neighbouring element. If the two filestreams are closed, output is suppressed.
Definition: binary_tree.cc:666
bool is_doc_enabled() const
Are we documenting?
Information for documentation of results: Directory and file number to enable output in the form RESL...
void find_neighbours()
Construct the neighbour lookup scheme.
Definition: binary_tree.cc:455
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
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: binary_tree.h:165
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
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
TreeRoot * Root_pt
Pointer to the root of the tree.
Definition: tree.h:270
BinaryTree * gteq_edge_neighbour(const int &direction, Vector< double > &s_in_neighbour, 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: binary_tree.cc:179
unsigned & number()
Number used (e.g.) for labeling output files.
static void setup_static_data()
Set up the static data, reflection schemes, etc.
Definition: binary_tree.cc:90
static DenseMatrix< int > Reflect
Reflection scheme: Reflect(direction,segment): Get mirror of segment in specified direction...
Definition: binary_tree.h:222
static Vector< std::string > Colour
Colours for neighbours in various directions.
Definition: binary_tree.h:207
int Level
Level of the Tree (level 0 = root)
Definition: tree.h:281
static DenseMatrix< bool > Is_adjacent
Array of direction/segment adjacency scheme: Is_adjacent(i_vertex,j_segment): Is vertex adjacent to s...
Definition: binary_tree.h:218
static char t char * s
Definition: cfortran.h:572
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: binary_tree.cc:547
unsigned self_test()
Self-test: Check all neighbours. Return success (0) if the maximum distance between corresponding poi...
Definition: binary_tree.cc:624
BinaryTreeForest()
Default constructor (empty and broken)
Definition: binary_tree.h:299
static Vector< int > Reflect_edge
Get opposite edge, e.g. Reflect_edge[L]=R.
Definition: binary_tree.h:214
static Vector< double > S_base
S_base(direction): Initial value for coordinate s on the edge indicated by direction (L/R) ...
Definition: binary_tree.h:211
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
void stick_all_tree_nodes_into_vector(Vector< Tree * > &)
Traverse and stick pointers to all "nodes" into Vector.
Definition: tree.cc:276
long number() const
Element number (for debugging/plotting)
unsigned self_test()
Self-test: Check all neighbours. Return success (0) if the maximum distance between corresponding poi...
Definition: binary_tree.cc:394
int Son_type
Son type (e.g. SW/SE/NW/NE in a quadtree)
Definition: tree.h:284
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:143
void resize(const unsigned long &n)
Definition: matrices.h:511