octree.h
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 #ifndef OOMPH_OCTREE_HEADER
31 #define OOMPH_OCTREE_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 //OOMPH-LIB headers
39 #include "tree.h"
40 #include "matrices.h"
41 
42 namespace oomph
43 {
44 
45 
46 
47 //====================================================================
48 // Namespace for OcTree directions
49 //====================================================================
50 namespace OcTreeNames
51 {
52 
53  /// \short Directions. OMEGA is used if a direction is undefined
54  /// in a certain context
55  enum{
56  LDB,RDB,LUB,RUB, // back octants/vertices
57  LDF,RDF,LUF,RUF, // front octants/vertices
58  LB,RB,DB,UB, // back edges
59  LD,RD,LU,RU, // side edges
60  LF,RF,DF,UF, // front edges
61  L,R,D,U,B,F, // faces
62  OMEGA=26};
63 }
64 
65 
66 // Forward definitions
67 class OcTreeRoot;
68 
69 //================================================================
70 /// OcTree class: Recursively defined, generalised octree.
71 ///
72 /// An OcTree has:
73 /// - a pointer to the object (of type RefineableQElement<3>) it represents
74 /// - Vector of pointers to its eight (LDB,RDB,...,RUF) sons (which are
75 /// octrees themselves). If the Vector of pointers to the sons has zero
76 /// length, the OcTree is a "leaf node" in the overall octree.
77 /// - a pointer to its father. If this pointer is NULL, the OcTree is the
78 /// the root node of the overall octree.
79 /// This data is stored in the Tree base class.
80 ///
81 /// The tree can also be part of a forest. If that is the case, the root
82 /// will have pointers to the roots of neighbouring octrees.
83 ///
84 /// The objects contained in the octree are assumed to be
85 /// (topologically) cubic elements whose geometry is
86 /// parametrised by local coordinates \f$ {\bf s} \in [-1,1]^3 \f$.
87 ///
88 /// The tree can be traversed while actions are being performed
89 /// at all of its "nodes" or only at the leaf "nodes".
90 ///
91 /// Finally, the leaf "nodes" can be split depending on criteria
92 /// defined by the object.
93 ///
94 /// Note that OcTrees are only generated by splitting existing
95 /// OcTrees. Therefore, the constructors are protected. The
96 /// only OcTree that "Joe User" can create is
97 /// the (derived) class OcTreeRoot.
98 //=================================================================
99 class OcTree : public virtual Tree
100 {
101 
102  public:
103 
104  /// \short Destructor. Note: Deleting a octree also deletes the
105  /// objects associated with all non-leaf nodes!
106  virtual ~OcTree() {}
107 
108 
109  /// Broken copy constructor
110  OcTree(const OcTree& dummy)
111  {
112  BrokenCopy::broken_copy("OcTree");
113  }
114 
115  /// Broken assignment operator
116  void operator=(const OcTree&)
117  {
118  BrokenCopy::broken_assign("OcTree");
119  }
120 
121 
122  /// \short Overload the function construct_son to ensure that the son
123  /// is a specific OcTree and not a general Tree.
125  Tree* const &father_pt, const int &son_type)
126  {
127  OcTree* temp_oc_pt = new OcTree(object_pt,father_pt,son_type);
128  return temp_oc_pt;
129  }
130 
131 
132 /// \short Find (pointer to) `greater-or-equal-sized face neighbour' in
133 /// given direction (L/R/U/D/F/B).
134 /// Another way of interpreting this is that we're looking for
135 /// the neighbour across the present element's face 'direction'.
136 /// The various arguments return additional information about the
137 /// size and relative orientation of the neighbouring octree.
138 /// To interpret these we use the following
139 /// <B>General convention:</B>
140 /// - Each face of the element that is represented by the octree
141 /// is parametrised by two (of the three)
142 /// local coordinates that parametrise the entire 3D element. E.g.
143 /// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
144 /// is parametrised by (s[0],s[2]); etc. We always identify the
145 /// in-face coordinate with the lower (3D) index with the subscript
146 /// _lo and the one with the larger (3D) index with the subscript _hi.
147 /// .
148 /// With this convention, the interpretation of the arguments is
149 /// as follows:
150 /// - The vector \c translate_s turns the index of the local coordinate
151 /// in the present octree into that of the neighbour. If there are no
152 /// rotations then \c translate_s[i] = i.
153 /// - In the present octree, the "south west" vertex of the face
154 /// between the present octree and its neighbour is located at
155 /// S_lo=-1, S_hi=-1. This point is located at the (3D) local
156 /// coordinates (\c s_sw[0], \c s_sw[1], \c s_sw[2])
157 /// in the neighbouring octree.
158 /// - ditto with s_ne: In the present octree, the "north east" vertex
159 /// of the face between the present octree and its neighbour is located at
160 /// S_lo=+1, S_hi=+1. This point is located
161 /// at the (3D) local coordinates (\c s_ne[0], \c s_ne[1], \c s_ne[2])
162 /// in the neighbouring octree.
163 /// - We're looking for a neighbour in the specified \c direction. When
164 /// viewed from the neighbouring octree, the face that separates
165 /// the present octree from its neighbour is the neighbour's face
166 /// \c face. If there's no rotation between the two octrees, this is a
167 /// simple reflection: For instance, if we're looking
168 /// for a neighhbour in the \c R [ight] \c direction, \c face will
169 /// be \c L [eft]
170 /// - \c diff_level <= 0 indicates the difference in refinement levels between
171 /// the two neighbours. If \c diff_level==0, the neighbour has the
172 /// same size as the current octree.
173  OcTree* gteq_face_neighbour(const int& direction,
174  Vector<unsigned> &translate_s,
175  Vector<double>& s_sw,
176  Vector<double>& s_ne,
177  int& face,
178  int& diff_level) const;
179 
180 
181 
182 /// \short Find (pointer to) `greater-or-equal-sized true edge neighbour' in
183 /// the given direction (LB,RB,DB,UB [the back edges],
184 /// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]).
185 ///
186 /// Another way of interpreting this is that we're looking for
187 /// the neighbour across the present element's edge 'direction'.
188 /// The various arguments return additional information about the
189 /// size and relative orientation of the neighbouring octree.
190 /// Each edge of the element that is represented by the octree
191 /// is parametrised by one (of the three) local coordinates that
192 /// parametrise the entire 3D element. E.g. the L[eft]B[ack] edge
193 /// is parametrised by s[1]; the "low" vertex of this edge
194 /// (located at the low value of this coordinate, i.e. at s[1]=-1)
195 /// is L[eft]D[own]B[ack]. The "high" vertex of this edge (located
196 /// at the high value of this coordinate, i.e. at s[1]=1) is
197 /// L[eft]U[p]B[ack]; etc
198 ///
199 /// The interpretation of the arguments is as follows:
200 /// - In a forest, an OcTree can have multiple edge neighbours
201 /// (across an edge where multiple trees meet). \c i_root_edge_neighbour
202 /// specifies which of these is used. Use this as "reverse communication":
203 /// First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
204 /// initialised to anything you want (zero, ideally). On return from
205 /// the fct, \c n_root_edge_neighour contains the total number of true
206 /// edge neighbours, so additional calls to the fct with
207 /// \c i_root_edge_neighbour>0 can be made until they've all been visited.
208 /// - The vector \c translate_s turns the index of the local coordinate
209 /// in the present octree into that of the neighbour. If there are no
210 /// rotations then \c translate_s[i] = i.
211 /// - The "low" vertex of the edge in the present octree
212 /// coincides with a certain vertex in the edge neighbour.
213 /// In terms of the neighbour's local coordinates, this point is
214 /// located at the (3D) local coordinates (\c s_lo[0], \c s_lo[1],
215 /// \c s_lo[2])
216 /// - ditto with s_hi: The "high" vertex of the edge in the present octree
217 /// coincides with a certain vertex in the edge neighbour.
218 /// In terms of the neighbour's local coordinates, this point is
219 /// located at the (3D) local coordinates (\c s_hi[0], \c s_hi[1],
220 /// \c s_hi[2])
221 /// - We're looking for a neighbour in the specified \c direction. When
222 /// viewed from the neighbouring octree, the edge that separates
223 /// the present octree from its neighbour is the neighbour's edge
224 /// \c edge. If there's no rotation between the two octrees, this is a
225 /// simple reflection: For instance, if we're looking
226 /// for a neighhbour in the \c DB \c direction, \c edge will
227 /// be \c UF.
228 /// - \c diff_level <= 0 indicates the difference in refinement levels between
229 /// the two neighbours. If \c diff_level==0, the neighbour has the
230 /// same size as the current octree.
231 /// .
232 /// \b Important: We're only looking for \b true edge neighbours
233 /// i.e. edge neigbours that are not also face neighbours. This is an
234 /// important difference to Samet's terminology. If the neighbour
235 /// in a certain direction is not a true edge neighbour, or if there
236 /// is no neighbour, then this function returns NULL.
237  OcTree* gteq_true_edge_neighbour(const int& direction,
238  const unsigned& i_root_edge_neighbour,
239  unsigned& nroot_edge_neighbour,
240  Vector<unsigned> &translate_s,
241  Vector<double>& s_lo,
242  Vector<double>& s_hi,
243  int& edge,
244  int& diff_level) const;
245 
246 
247  /// \short Self-test: Check all neighbours. Return success (0)
248  /// if the max. distance between corresponding points in the
249  /// neighbours is less than the tolerance specified in the
250  /// static value Tree::Max_neighbour_finding_tolerance.
251  unsigned self_test();
252 
253  /// \short Setup the static data, rotation and reflection schemes, etc
254  static void setup_static_data();
255 
256 
257  /// \short Doc/check all face neighbours of octree (nodes) contained in the
258  /// Vector forest_node_pt. Output into neighbours_file which can
259  /// be viewed from tecplot with OcTreeNeighbours.mcr
260  /// Neighbour info and errors are displayed on
261  /// neighbours_txt_file. Finally, compute the max. error between
262  /// vertices when viewed from neighhbouring element.
263  /// If the two filestreams are closed, output is suppressed.
264  static void doc_face_neighbours(Vector<Tree*> forest_nodes_pt,
265  std::ofstream& neighbours_file,
266  std::ofstream& neighbours_txt_file,
267  double& max_error);
268 
269 
270  /// \short Doc/check all true edge neighbours of octree (nodes) contained
271  /// in the Vector forest_node_pt. Output into neighbours_file which can
272  /// be viewed from tecplot with OcTreeNeighbours.mcr
273  /// Neighbour info and errors are displayed on
274  /// neighbours_txt_file. Finally, compute the max. error between
275  /// vertices when viewed from neighhbouring element.
276  /// If the two filestreams are closed, output is suppressed.
277  static void doc_true_edge_neighbours(Vector<Tree*> forest_nodes_pt,
278  std::ofstream& neighbours_file,
279  std::ofstream& no_true_edge_file,
280  std::ofstream& neighbours_txt_file,
281  double& max_error);
282 
283  /// \short If an edge is bordered by the nodes whose local numbers
284  /// are n1 and n2 in an element with nnode1d nodes along each coordinate
285  /// direction, then this edge is shared by two faces. This function
286  /// takes one of these faces as the argument \c face and returns the
287  /// other one. (\c face is a direction in the set U,D,F,B,L,R).
288  static int get_the_other_face(const unsigned& n1, const unsigned& n2,
289  const unsigned& nnode1d, const int& face);
290 
291  /// \short Return the local node number of given vertex
292  /// [LDB,RDB,...] in an element with nnode1d nodes in each
293  /// coordinate direction
294  static unsigned vertex_to_node_number(const int& vertex,
295  const unsigned& nnode1d);
296 
297 
298  /// \short Return the vertex [LDB,RDB,...] of local (vertex) node n
299  /// in an element with nnode1d nodes in each coordinate direction.
300  static int node_number_to_vertex(const unsigned& n,
301  const unsigned& nnode1d);
302 
303 
304  /// \short If U[p] becomes new_up and R[ight] becomes new_right then the
305  /// direction vector \c dir becomes rotate(new_up, new_right, dir)
306  static Vector<int> rotate(const int& new_up, const int& new_right,
307  const Vector<int>& dir);
308 
309 
310  /// \short If U[p] becomes new_up and R[ight] becomes new_right
311  /// then the direction \c dir becomes \c rotate(new_up, new_right, dir)
312  static int rotate(const int& new_up, const int& new_right,
313  const int& dir);
314 
315 
316 
317  /// Translate (enumerated) directions into strings
319 
320  /// Get opposite face, e.g. Reflect_face[L]=R
322 
323  /// Get opposite edge, e.g. Reflect_edge[DB]=UF
325 
326  /// Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF
328 
329  /// \short \c Vector of vectors containing the two vertices for each edge,
330  /// e.g. \c Vertex_at_end_of_edge[LU][0]=LUB and
331  /// \c Vertex_at_end_of_edge[LU][1]=LUF.
333 
334  /// \short Each vector representing a direction can be translated into
335  /// a direction, either a son type (vertex), a face or an edge.
336  /// E.g. : Vector_to_direction[(1,-1,1)]=RDF, Vector_to_direction[(0,1,0)]=U
337  static std::map<Vector<int>, int > Vector_to_direction;
338 
339  /// \short For each direction, i.e. a son_type (vertex), a face or an
340  /// edge, this defines a vector that indicates this direction.
341  /// E.g : Direction_to_vector[RDB]=(1,-1,-1), Direction_to_vector[U]=(0,1,0)
343 
344 
345  /// \short Storage for the up/right-equivalents corresponding to two
346  /// pairs of vertices along an element edge:
347  /// - The first pair contains
348  /// -# the vertex in the reference element
349  /// -# the corresponding vertex in the edge neighbour (i.e. the
350  /// vertex in the edge neighbour that is located at the same
351  /// position as that first vertex).
352  /// .
353  /// - The second pair contains
354  /// -# the vertex at the other end of the edge in the reference element
355  /// -# the corresponding vertex in the edge neighbour.
356  /// .
357  /// .
358  /// These two pairs completely define the relative rotation
359  /// between the reference element and its edge neighbour. The map
360  /// returns a pair which contains the up_equivalent and the
361  /// right_equivalent of the edge neighbour, i.e. it tells us
362  /// which direction in the edge neighbour coincides with the
363  /// up (or right) direction in the reference element.
364  static std::map<std::pair<std::pair<int,int>,std::pair<int,int> >,
366 
367 
368  protected:
369 
370  /// Default constructor (empty and broken)
372  {
373  throw OomphLibError(
374  "Don't call empty constructor for OcTree!",
375  OOMPH_CURRENT_FUNCTION,
376  OOMPH_EXCEPTION_LOCATION);
377  }
378 
379  /// \short Constructor for empty (root) tree:
380  /// no father, no sons; just pass a pointer to its object
381  /// (a RefineableQElement<3>). This is
382  /// protected because OcTrees can only be created internally,
383  /// during the split operation. Only OcTreeRoots can be
384  /// created externally.
385  OcTree(RefineableElement* const &object_pt) : Tree(object_pt) {}
386 
387 
388  /// \short Constructor for tree that has a father: Pass it the pointer
389  /// to its object, the pointer to its father and tell it what type
390  /// of son (LDB,RDB,...) it is.
391  /// Protected because OcTrees can only be created internally,
392  /// during the split operation. Only OcTreeRoots can be
393  /// created externally.
395  Tree* const &father_pt, const int& son_type) :
396  Tree(object_pt,father_pt,son_type) {}
397 
398 
399  /// Bool indicating that static member data has been setup
401 
402 
403  private:
404 
405  /// \short Find `greater-or-equal-sized face neighbour' in given direction
406  /// (L/R/U/D/B/F).
407  ///
408  /// This is an auxiliary routine which allows neighbour finding in adjacent
409  /// octrees. Needs to keep track of the maximum level to which
410  /// search is performed because in the presence of OcTree forests,
411  /// the search isn't purely recursive.
412  ///
413  /// Parameters:
414  /// - direction: (L/R/U/D/B/F) Direction in which neighbour has to be found.
415  /// - s_difflo/s_diffhi: Offset of left/down/back vertex from
416  /// corresponding vertex in
417  /// neighbour. Note that this is input/output as it needs to be incremented/
418  /// decremented during the recursive calls to this function.
419  /// - diff_level <= 0 indicates the difference in octree levels
420  /// between the current element and its neighbour.
421  /// - max_level is the maximum level to which the neighbour search is
422  /// allowed to proceed. This is necessary because in a forest,
423  /// the neighbour search isn't based on pure recursion.
424  /// - orig_root_pt identifies the root node of the element whose
425  /// neighbour we're really trying to find by all these recursive calls.
426  OcTree* gteq_face_neighbour(const int& direction,
427  double& s_difflo,
428  double& s_diffhi,
429  int& diff_level ,
430  int max_level,
431  OcTreeRoot* orig_root_pt) const;
432 
433 
434  /// \short Find `greater-or-equal-sized edge neighbour' in given direction
435  /// (LB,RB,DB,UB [the back edges],
436  /// LD,RD,LU,RU [the side edges], LF,RF,DF,UF [the front edges]).
437  ///
438  /// This is an auxiliary routine which allows neighbour finding in adjacent
439  /// octrees. Needs to keep track of the maximum level to which
440  /// search is performed because in the presence of OcTree forests,
441  /// the search isn't purely recursive.
442  ///
443  /// Parameters:
444  /// - direction: (LB/RB/...) Direction in which neighbour has to be found.
445  /// - In a forest, an OcTree can have multiple edge neighbours
446  /// (across an edge where multiple trees meet). \c i_root_edge_neighbour
447  /// specifies which of these is used. Use this as "reverse communication":
448  /// First call with \c i_root_edge_neighbour=0 and \c n_root_edge_neighour
449  /// initialised to anything you want (zero, ideally). On return from
450  /// the fct, \c n_root_edge_neighour contains the total number of true
451  /// edge neighbours, so additional calls to the fct with
452  /// \c i_root_edge_neighbour>0 can be made until they've all been visited.
453  /// - s_diff: Offset of the edge's "low" vertex from
454  /// corresponding vertex in
455  /// neighbour. Note that this is input/output as it needs to be incremented/
456  /// decremented during the recursive calls to this function.
457  /// - diff_level <= 0 indicates the difference in octree levels
458  /// between the current element and its neighbour.
459  /// - max_level is the maximum level to which the neighbour search is
460  /// allowed to proceed. This is necessary because in a forest,
461  /// the neighbour search isn't based on pure recursion.
462  /// - orig_root_pt identifies the root node of the element whose
463  /// neighbour we're really trying to find by all these recursive calls.
464  /// .
465  /// \b Note: some of the auxiliary information may be incorrect if
466  /// the neighbour is not a true edge neighbour. We don't care because
467  /// we're not dealing with those!
468  OcTree* gteq_edge_neighbour(const int& direction,
469  const unsigned& i_root_edge_neighbour,
470  unsigned& nroot_edge_neighbour,
471  double& s_diff,
472  int& diff_level ,
473  int max_level,
474  OcTreeRoot* orig_root_pt) const;
475 
476 
477  /// \short Is the edge neighbour (for edge "edge") specified via the pointer
478  /// also a face neighbour for one of the two adjacent faces?
479  bool edge_neighbour_is_face_neighbour(const int& edge,
480  OcTree* edge_neighb_pt) const;
481 
482 
483 
484 
485  /// \short This constructs the rotation matrix of the rotation around the
486  /// axis \c axis with an angle of \c angle*90
487  static void construct_rotation_matrix(int& axis, int& angle,
488  DenseMatrix<int>& mat);
489 
490  /// Helper function: Performs the operation : vect2 = mat*vect1
491  static void mult_mat_vect(const DenseMatrix<int>& mat,
492  const Vector<int>& vect1,
493  Vector<int>& vect2);
494 
495  /// Helper function: Performs the operation : mat3=mat1*mat2
496  static void mult_mat_mat(const DenseMatrix<int>& mat1,
497  const DenseMatrix<int>& mat2,
498  DenseMatrix<int>& mat3);
499 
500 
501  /// \short Returns the vector of the coordinate directions
502  /// of vertex node number n in an element with nnode1d element per
503  /// dimension.
504  static Vector<int> vertex_node_to_vector(const unsigned& n,
505  const unsigned& nnode1d);
506 
507 
508 
509  /// Entry in rotation matrix: cos(i*90)
511 
512  /// Entry in rotation matrix sin(i*90)
514 
515 
516  /// \short Array of direction/octant adjacency scheme:
517  /// Is_adjacent(direction,octant): Is face/edge \c direction
518  /// adjacent to octant \c octant ? (Table in Samet's book)
520 
521  /// \short Reflection scheme: Reflect(direction,octant): Get mirror
522  /// of octant/edge in specified direction. E.g. Reflect(LDF,L)=RDF
524 
525  /// \short Determine common face of edges or octants.
526  /// Slightly bizarre lookup scheme from Samet's book.
528 
529  /// Colours for neighbours in various directions
531 
532  /// \short s_base(i,direction): Initial value for coordinate s[i] on
533  /// the face indicated by direction (L/R/U/D/F/B)
535 
536  /// \short Each face of the RefineableQElement<3> that is represented
537  /// by the octree is parametrised by two (of the three)
538  /// local coordinates that parametrise the entire 3D element. E.g.
539  /// the B[ack] face is parametrised by (s[0], s[1]); the D[own] face
540  /// is parametrised by (s[0],s[2]); etc. We always identify the
541  /// in-face coordinate with the lower (3D) index with the subscript
542  /// \c _lo and the one with the larger (3D) index with the subscript \c _hi.
543  /// Here we set up the translation scheme between the 2D in-face
544  /// coordinates (s_lo,s_hi) and the corresponding 3D coordinates:
545  /// If we're located on face \c face [L/R/F/B/U/D], then
546  /// an increase in s_lo from -1 to +1 corresponds to a change
547  /// of \c s_steplo(i,face) in the 3D coordinate \c s[i].
549 
550  /// \short If we're located on face \c face [L/R/F/B/U/D], then
551  /// an increase in s_hi from -1 to +1 corresponds to a change
552  /// of \c s_stephi(i,face) in the 3D coordinate \ s[i].
553  /// [Read the discussion of \c s_steplo for an explanation of
554  /// the subscripts \c _hi and \c _lo.]
556 
557  /// \short Relative to the left/down/back vertex in any (father) octree, the
558  /// corresponding vertex in the son specified by \c son_octant has an offset.
559  /// If we project the son_octant's left/down/back vertex onto the
560  /// father's face \c face, it is located at the in-face coordinate
561  /// \c s_lo = h/2 \c S_directlo(face,son_octant). [See discussion of
562  /// \c s_steplo for an explanation of the subscripts \c _hi and \c _lo.]
564 
565  /// \short Relative to the left/down/back vertex in any (father) octree, the
566  /// corresponding vertex in the son specified by \c son_octant has an offset.
567  /// If we project the son_octant's left/down/back vertex onto the
568  /// father's face \c face, it is located at the in-face coordinate
569  /// \c s_hi = h/2 \c S_directlhi(face,son_octant). [See discussion of
570  /// \c s_steplo for an explanation of the subscripts \c _hi and \c _lo.]
572 
573  /// \short S_base_edge(i,edge): Initial value for coordinate s[i] on
574  /// the specified edge (LF/RF/...).
576 
577  /// \short Each edge of the RefineableQElement<3> that is represented
578  /// by the octree is parametrised by one (of the three)
579  /// local coordinates that parametrise the entire 3D element.
580  /// If we're located on edge \c edge [DB,UB,...], then
581  /// an increase in s from -1 to +1 corresponds to a change
582  /// of \c s_step_edge(i,edge) in the 3D coordinates \c s[i].
584 
585  /// \short Relative to the left/down/back vertex in any (father) octree, the
586  /// corresponding vertex in the son specified by \c son_octant has an offset.
587  /// If we project the son_octant's left/down/back vertex onto the
588  /// father's edge \c edge, it is located at the in-face coordinate
589  /// \c s_lo = h/2 \c S_direct_edge(edge,son_octant).
591 
592 
593 };
594 
595 
596 
597 //===================================================================
598 /// OcTreeRoot is a OcTree that forms the root of a (recursive)
599 /// octree. The "root node" is special as it holds additional
600 /// information about its neighbours and their relative
601 /// rotation (inside a OcTreeForest).
602 //==================================================================
603 class OcTreeRoot : public virtual OcTree, public virtual TreeRoot
604 {
605  private:
606 
607  /// \short Map of pointers to the edge-neighbouring [Oc]TreeRoots:
608  /// Edge_neighbour_pt[direction] is Vector to the pointers to the
609  /// [Oc]TreeRoot's edge neighbours in the (enumerated) (edge) direction.
610  std::map<int,Vector<TreeRoot*> > Edge_neighbour_pt;
611 
612  /// \short Map giving the Up equivalent of the neighbour specified by
613  /// pointer: When viewed from the current octree's neighbour,
614  /// our up direction is the neighbour's Up_equivalent[neighbour_pt]
615  /// direction. If there's no rotation, this map contains the identify
616  /// so that, e.g. \c Up_equivalent[neighbour_pt]=U (read as: "in my
617  /// neighbour, my Up is its Up"). If the neighbour is rotated
618  /// by 180 degrees relative to the current octree(around the back-front axis),
619  /// say, we have \c Up_equivalent[neighbour_pt]=D (read as: "in my
620  /// neighbour, my Up is its Down"); etc.
621  std::map<TreeRoot*,int> Up_equivalent;
622 
623  /// \short Map giving the Right equivalent of the neighbour specified by
624  /// pointer: When viewed from the current octree's neighbour,
625  /// our right direction is the neighbour's right_equivalent[neighbour_pt]
626  /// direction. If there's no rotation, this map contains the identify
627  /// so that, e.g. \c Right_equivalent[neighbour_pt]=R (read as: "in my
628  /// neighbour, my Right is its Right").
629  std::map<TreeRoot*,int> Right_equivalent;
630 
631 
632  public:
633 
634  /// \short Constructor for the root octree: Pass pointer to the
635  /// RefineableQElement<3> that is represented by the OcTree.
637  OcTree(object_pt), TreeRoot(object_pt)
638  {
639 
640 #ifdef PARANOID
641  // Check that static member data has been set up
643  {
644  std::string error_message =
645  "Static member data hasn't been setup yet.\n";
646  error_message +=
647  "Call OcTree::setup_static_data() before creating\n";
648  error_message += "any OcTreeRoots\n";
649 
650  throw OomphLibError(error_message,
651  OOMPH_CURRENT_FUNCTION,
652  OOMPH_EXCEPTION_LOCATION);
653  }
654 #endif
655 
656  }
657 
658 
659  /// Broken copy constructor
660  OcTreeRoot(const OcTreeRoot& dummy) : TreeRoot(dummy)
661  {
662  BrokenCopy::broken_copy("OcTreeRoot");
663  }
664 
665  /// Broken assignment operator
666  void operator=(const OcTreeRoot&)
667  {
668  BrokenCopy::broken_assign("OcTreeRoot");
669  }
670 
671 
672  /// \short Return vector of pointers to the edge-neighbouring TreeRoots
673  /// in the (enumerated) (edge) direction.
674  Vector<TreeRoot*> edge_neighbour_pt(const unsigned& edge_direction)
675  {
676 
677 #ifdef PARANOID
678  using namespace OcTreeNames;
679  if ((edge_direction!=LB)&&(edge_direction!=RB)&&(edge_direction!=DB)&&
680  (edge_direction!=UB)&&
681  (edge_direction!=LD)&&(edge_direction!=RD)&&(edge_direction!=LU)&&
682  (edge_direction!=RU)&&
683  (edge_direction!=LF)&&(edge_direction!=RF)&&(edge_direction!=DF)&&
684  (edge_direction!=UF))
685  {
686  std::ostringstream error_stream;
687  error_stream << "Wrong edge_direction: " << Direct_string[edge_direction]
688  << std::endl;
689  throw OomphLibError(error_stream.str(),
690  OOMPH_CURRENT_FUNCTION,
691  OOMPH_EXCEPTION_LOCATION);
692  }
693 #endif
694 
695  return Edge_neighbour_pt[edge_direction];
696  }
697 
698 
699  /// \short Return number of edge-neighbouring OcTreeRoot
700  /// in the (enumerated) (edge) direction.
701  unsigned nedge_neighbour(const unsigned& edge_direction)
702  {
703 
704 #ifdef PARANOID
705  using namespace OcTreeNames;
706  if ((edge_direction!=LB)&&(edge_direction!=RB)&&(edge_direction!=DB)&&
707  (edge_direction!=UB)&&
708  (edge_direction!=LD)&&(edge_direction!=RD)&&(edge_direction!=LU)&&
709  (edge_direction!=RU)&&
710  (edge_direction!=LF)&&(edge_direction!=RF)&&(edge_direction!=DF)&&
711  (edge_direction!=UF))
712  {
713  std::ostringstream error_stream;
714  error_stream << "Wrong edge_direction: " << Direct_string[edge_direction]
715  << std::endl;
716  throw OomphLibError(error_stream.str(),
717  OOMPH_CURRENT_FUNCTION,
718  OOMPH_EXCEPTION_LOCATION);
719  }
720 #endif
721  return Edge_neighbour_pt[edge_direction].size();
722  }
723 
724  /// \short Add pointer to the edge-neighbouring [Oc]TreeRoot
725  /// in the (enumerated) (edge) direction -- maintains uniqueness
726  void add_edge_neighbour_pt(TreeRoot* oc_tree_root_pt,
727  const unsigned& edge_direction)
728  {
729 
730 #ifdef PARANOID
731  using namespace OcTreeNames;
732  if ((edge_direction!=LB)&&(edge_direction!=RB)&&(edge_direction!=DB)&&
733  (edge_direction!=UB)&&
734  (edge_direction!=LD)&&(edge_direction!=RD)&&(edge_direction!=LU)&&
735  (edge_direction!=RU)&&
736  (edge_direction!=LF)&&(edge_direction!=RF)&&(edge_direction!=DF)&&
737  (edge_direction!=UF))
738  {
739  std::ostringstream error_stream;
740  error_stream << "Wrong edge_direction: " << Direct_string[edge_direction]
741  << std::endl;
742  throw OomphLibError(error_stream.str(),
743  OOMPH_CURRENT_FUNCTION,
744  OOMPH_EXCEPTION_LOCATION);
745  }
746 #endif
747 
749  find(Edge_neighbour_pt[edge_direction].begin(),
750  Edge_neighbour_pt[edge_direction].end(),
751  oc_tree_root_pt);
752  if (it==Edge_neighbour_pt[edge_direction].end())
753  {
754  Edge_neighbour_pt[edge_direction].push_back(oc_tree_root_pt);
755  }
756  }
757 
758 
759  /// \short Return up equivalent of the neighbours specified by
760  /// pointer: When viewed from the current octree's neighbour,
761  /// our up direction is the neighbour's Up_equivalent[neighbour_pt]
762  /// direction. If there's no rotation, this map contains the identify
763  /// so that, e.g. \c Up_equivalent[neighbour_pt]=U (read as: "in my
764  /// neighbour, my Up is its Up"). If the neighbour is rotated
765  /// by 180 degrees relative to the current octree (around the back-front axis)
766  /// say, we have \c Up_equivalent[neighbour_pt]=D (read as: "in my
767  /// neighbour, my Up is its Down"); etc. Returns OMEGA if the Octree
768  /// specified by the pointer argument is not a neighbour.
769  int up_equivalent(TreeRoot* tree_root_pt)
770  {
771  if (direction_of_neighbour(tree_root_pt)==OMEGA)
772  {
773  return OMEGA;
774  }
775  else
776  {
777  return Up_equivalent[tree_root_pt];
778  }
779  }
780 
781 
782  /// \short Set up equivalent of the neighbours specified by
783  /// pointer: When viewed from the current octree's neighbour,
784  /// our up direction is the neighbour's Up_equivalent[neighbour_pt]
785  /// direction. If there's no rotation, this map contains the identify
786  /// so that, e.g. \c Up_equivalent[neighbour_pt]=U (read as: "in my
787  /// neighbour, my Up is its Up"). If the neighbour is rotated
788  /// by 180 degrees relative to the current octree (around the back-front axis)
789  /// say, we have \c Up_equivalent[neighbour_pt]=D (read as: "in my
790  /// neighbour, my Up is its Down"); etc.
791  void set_up_equivalent(TreeRoot* tree_root_pt, const int& dir)
792  {
793  Up_equivalent[tree_root_pt]=dir;
794  }
795 
796 
797  /// \short The same thing as up_equivalent, but for the right direction:
798  /// When viewed from the current octree neighbour, our
799  /// right direction is the neighbour's Right_equivalent[neighbour_pt]
800  /// direction. Returns OMEGA if the Octree specified by the pointer
801  /// argument is not a neighbour.
802  int right_equivalent(TreeRoot* tree_root_pt)
803  {
804  if (direction_of_neighbour(tree_root_pt)==OMEGA)
805  {
806  return OMEGA;
807  }
808  else
809  {
810  return Right_equivalent[tree_root_pt];
811  }
812  }
813 
814  /// \short The same thing as up_equivalent, but for the right direction:
815  /// When viewed from the current octree neighbour, our
816  /// right direction is the neighbour's Right_equivalent[neighbour_pt]
817  /// direction.
818  void set_right_equivalent(TreeRoot* tree_root_pt, const int& dir)
819  {
820  Right_equivalent[tree_root_pt]=dir;
821  }
822 
823  /// \short If octree_root_pt is a neighbour, return the direction
824  /// [faces L/R/F/B/U/D or edges DB/UP/...] in which it is found,
825  /// otherwise return OMEGA
826  int direction_of_neighbour(TreeRoot* octree_root_pt)
827  {
828  using namespace OcTreeNames;
829 
830  if(Neighbour_pt[U]==octree_root_pt) {return U;}
831  if(Neighbour_pt[D]==octree_root_pt) {return D;}
832  if(Neighbour_pt[L]==octree_root_pt) {return L;}
833  if(Neighbour_pt[R]==octree_root_pt) {return R;}
834  if(Neighbour_pt[F]==octree_root_pt) {return F;}
835  if(Neighbour_pt[B]==octree_root_pt) {return B;}
836 
837  if(Neighbour_pt[LB]==octree_root_pt) {return LB;}
838  if(Neighbour_pt[RB]==octree_root_pt) {return RB;}
839  if(Neighbour_pt[DB]==octree_root_pt) {return DB;}
840  if(Neighbour_pt[UB]==octree_root_pt) {return UB;}
841 
842  if(Neighbour_pt[LD]==octree_root_pt) {return LD;}
843  if(Neighbour_pt[RD]==octree_root_pt) {return RD;}
844  if(Neighbour_pt[LU]==octree_root_pt) {return LU;}
845  if(Neighbour_pt[RU]==octree_root_pt) {return RU;}
846 
847  if(Neighbour_pt[LF]==octree_root_pt) {return LF;}
848  if(Neighbour_pt[RF]==octree_root_pt) {return RF;}
849  if(Neighbour_pt[DF]==octree_root_pt) {return DF;}
850  if(Neighbour_pt[UF]==octree_root_pt) {return UF;}
851 
852 
853  // Search over all edge neighbours
854  for (int dir=LB;dir<=UF;dir++)
855  {
856  Vector<TreeRoot*> edge_neigh_pt=this->edge_neighbour_pt(dir);
857  unsigned n_neigh=edge_neigh_pt.size();
858  for (unsigned e=0;e<n_neigh;e++)
859  {
860  if (edge_neigh_pt[e]==octree_root_pt)
861  {
862  return dir;
863  }
864  }
865  }
866 
867 
868  //If we get here, it's not a neighbour
869  return OMEGA;
870  }
871 
872 };
873 
874 
875 
876 //////////////////////////////////////////////////////////////////////////
877 //////////////////////////////////////////////////////////////////////////
878 //////////////////////////////////////////////////////////////////////////
879 
880 
881 //================================================================
882 /// An OcTreeForest consists of a collection of OcTreeRoots.
883 /// Each member tree can have neighbours to its L/R/U/D/F/B and
884 /// DB/UP/... and the orientation of their compasses can differ,
885 /// allowing for complex, unstructured meshes.
886 //=================================================================
887 class OcTreeForest : public TreeForest
888 {
889 
890  public:
891 
892  /// \short Constructor for OcTree forest: Pass Vector of
893  /// (pointers to) trees.
894  OcTreeForest(Vector<TreeRoot* >& trees_pt);
895 
896  /// Default constructor (empty and broken)
898  {
899  throw OomphLibError("Don't call empty constructor for OcTreeForest!",
900  OOMPH_CURRENT_FUNCTION,
901  OOMPH_EXCEPTION_LOCATION);
902  }
903 
904  /// Broken copy constructor
905  OcTreeForest(const OcTreeForest& dummy)
906  {
907  BrokenCopy::broken_copy("OcTreeForest");
908  }
909 
910  /// Broken assignment operator
911  void operator=(const OcTreeForest&)
912  {
913  BrokenCopy::broken_assign("OcTreeForest");
914  }
915 
916  /// \short Destructor: Delete the constituent octrees (and thus
917  /// the associated objects!)
918  virtual ~OcTreeForest() {}
919 
920 
921 
922  /// \short Document and check all the neighbours of all the nodes
923  /// in the forest. DocInfo object specifies the output directory
924  /// and file numbers for the various files. If \c doc_info.disable_doc()
925  /// has been called, no output is created.
926  void check_all_neighbours(DocInfo &doc_info);
927 
928 /// \short Open output files that will store any hanging nodes in
929  /// the forest and return a vector of the streams.
930  void open_hanging_node_files(DocInfo &doc_info,
931  Vector<std::ofstream*> &output_stream);
932 
933  /// \short Self-test: Check all neighbours. Return success (0)
934  /// if the max. distance between corresponding points in the
935  /// neighbours is less than the tolerance specified in the
936  /// static value Tree::Max_neighbour_finding_tolerance.
937  unsigned self_test();
938 
939  /// \short Return pointer to i-th OcTree in forest
940  /// (Performs a dynamic cast from the TreeRoot to a
941  /// OcTreeRoot).
942  OcTreeRoot* octree_pt(const unsigned &i) const
943  {return dynamic_cast<OcTreeRoot*>(Trees_pt[i]);}
944 
945  /// \short Given the number i of the root octree in this forest, return
946  /// pointer to its face neighbour in the specified direction. NULL
947  /// if neighbour doesn't exist. (This does the dynamic cast
948  /// from a TreeRoot to a OcTreeRoot internally).
949  OcTreeRoot* oc_face_neigh_pt(const unsigned &i, const int &direction)
950  {
951 #ifdef PARANOID
952  using namespace OcTreeNames;
953  if ((direction!=U)&&(direction!=D)&&
954  (direction!=F)&&(direction!=B)&&
955  (direction!=L)&&(direction!=R))
956  {
957  std::ostringstream error_stream;
958  error_stream << "Wrong edge_direction: "
959  << OcTree::Direct_string[direction]
960  << std::endl;
961  throw OomphLibError(error_stream.str(),
962  OOMPH_CURRENT_FUNCTION,
963  OOMPH_EXCEPTION_LOCATION);
964  }
965 #endif
966  return dynamic_cast<OcTreeRoot*>(Trees_pt[i]->neighbour_pt(direction));
967  }
968 
969  /// \short Given the number i of the root octree in this forest, return
970  /// the vector of pointers to the true edge neighbours in the specified
971  /// (edge) direction.
972  Vector<TreeRoot*> oc_edge_neigh_pt(const unsigned &i, const int &direction)
973  {
974  // Note: paranoia check is done in edge_neighbour_pt
975  return dynamic_cast<OcTreeRoot*>(
976  Trees_pt[i])->edge_neighbour_pt(direction);
977  }
978 
979 
980 private:
981 
982 
983  /// Construct the rotation schemes
984  void construct_up_right_equivalents();
985 
986  /// Construct the neighbour scheme
987  void find_neighbours();
988 
989 
990 };
991 
992 }
993 
994 #endif
995 
static std::map< std::pair< std::pair< int, int >, std::pair< int, int > >, std::pair< int, int > > Up_and_right_equivalent_for_pairs_of_vertices
Storage for the up/right-equivalents corresponding to two pairs of vertices along an element edge: ...
Definition: octree.h:365
static Vector< int > Reflect_face
Get opposite face, e.g. Reflect_face[L]=R.
Definition: octree.h:321
void operator=(const OcTreeRoot &)
Broken assignment operator.
Definition: octree.h:666
void set_right_equivalent(TreeRoot *tree_root_pt, const int &dir)
The same thing as up_equivalent, but for the right direction: When viewed from the current octree nei...
Definition: octree.h:818
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
static void setup_static_data()
Setup the static data, rotation and reflection schemes, etc.
Definition: octree.cc:866
OcTree * gteq_face_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_sw, Vector< double > &s_ne, int &face, int &diff_level) const
Find (pointer to) `greater-or-equal-sized face neighbour' in given direction (L/R/U/D/F/B). Another way of interpreting this is that we're looking for the neighbour across the present element's face 'direction'. The various arguments return additional information about the size and relative orientation of the neighbouring octree. To interpret these we use the following General convention:
Definition: octree.cc:3167
static unsigned vertex_to_node_number(const int &vertex, const unsigned &nnode1d)
Return the local node number of given vertex [LDB,RDB,...] in an element with nnode1d nodes in each c...
Definition: octree.cc:289
virtual ~OcTree()
Destructor. Note: Deleting a octree also deletes the objects associated with all non-leaf nodes! ...
Definition: octree.h:106
static std::map< Vector< int >, int > Vector_to_direction
Each vector representing a direction can be translated into a direction, either a son type (vertex)...
Definition: octree.h:337
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:102
Tree * construct_son(RefineableElement *const &object_pt, Tree *const &father_pt, const int &son_type)
Overload the function construct_son to ensure that the son is a specific OcTree and not a general Tre...
Definition: octree.h:124
static Vector< int > vertex_node_to_vector(const unsigned &n, const unsigned &nnode1d)
Returns the vector of the coordinate directions of vertex node number n in an element with nnode1d el...
Definition: octree.cc:245
static DenseMatrix< double > S_base_edge
S_base_edge(i,edge): Initial value for coordinate s[i] on the specified edge (LF/RF/...).
Definition: octree.h:575
static Vector< std::string > Colour
Colours for neighbours in various directions.
Definition: octree.h:530
Information for documentation of results: Directory and file number to enable output in the form RESL...
cstr elem_len * i
Definition: cfortran.h:607
int up_equivalent(TreeRoot *tree_root_pt)
Return up equivalent of the neighbours specified by pointer: When viewed from the current octree's ne...
Definition: octree.h:769
OcTree * gteq_true_edge_neighbour(const int &direction, const unsigned &i_root_edge_neighbour, unsigned &nroot_edge_neighbour, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level) const
Find (pointer to) `greater-or-equal-sized true edge neighbour' in the given direction (LB...
Definition: octree.cc:3416
OcTree * gteq_edge_neighbour(const int &direction, const unsigned &i_root_edge_neighbour, unsigned &nroot_edge_neighbour, double &s_diff, int &diff_level, int max_level, OcTreeRoot *orig_root_pt) const
Find `greater-or-equal-sized edge neighbour' in given direction (LB,RB,DB,UB [the back edges]...
Definition: octree.cc:3832
static DenseMatrix< double > S_steplo
Each face of the RefineableQElement<3> that is represented by the octree is parametrised by two (of the...
Definition: octree.h:548
bool edge_neighbour_is_face_neighbour(const int &edge, OcTree *edge_neighb_pt) const
Is the edge neighbour (for edge "edge") specified via the pointer also a face neighbour for one of th...
Definition: octree.cc:2591
std::map< int, Vector< TreeRoot * > > Edge_neighbour_pt
Map of pointers to the edge-neighbouring [Oc]TreeRoots: Edge_neighbour_pt[direction] is Vector to the...
Definition: octree.h:610
e
Definition: cfortran.h:575
static Vector< Vector< int > > Vertex_at_end_of_edge
Vector of vectors containing the two vertices for each edge, e.g. Vertex_at_end_of_edge[LU][0]=LUB an...
Definition: octree.h:332
unsigned nedge_neighbour(const unsigned &edge_direction)
Return number of edge-neighbouring OcTreeRoot in the (enumerated) (edge) direction.
Definition: octree.h:701
static DenseMatrix< double > S_base
s_base(i,direction): Initial value for coordinate s[i] on the face indicated by direction (L/R/U/D/F/...
Definition: octree.h:534
static Vector< int > rotate(const int &new_up, const int &new_right, const Vector< int > &dir)
If U[p] becomes new_up and R[ight] becomes new_right then the direction vector dir becomes rotate(new...
Definition: octree.cc:604
static Vector< int > Cosi
Entry in rotation matrix: cos(i*90)
Definition: octree.h:510
static Vector< int > Reflect_vertex
Get opposite vertex, e.g. Reflect_vertex[LDB]=RUF.
Definition: octree.h:327
double angle(const Vector< double > &a, const Vector< double > &b)
Get the angle between two vector.
Definition: Vector.h:307
static bool Static_data_has_been_setup
Bool indicating that static member data has been setup.
Definition: octree.h:400
static DenseMatrix< double > S_stephi
If we're located on face face [L/R/F/B/U/D], then an increase in s_hi from -1 to +1 corresponds to a ...
Definition: octree.h:555
static Vector< int > Reflect_edge
Get opposite edge, e.g. Reflect_edge[DB]=UF.
Definition: octree.h:324
static DenseMatrix< int > Reflect
Reflection scheme: Reflect(direction,octant): Get mirror of octant/edge in specified direction...
Definition: octree.h:523
OcTreeRoot(const OcTreeRoot &dummy)
Broken copy constructor.
Definition: octree.h:660
static void doc_true_edge_neighbours(Vector< Tree * > forest_nodes_pt, std::ofstream &neighbours_file, std::ofstream &no_true_edge_file, std::ofstream &neighbours_txt_file, double &max_error)
Doc/check all true edge neighbours of octree (nodes) contained in the Vector forest_node_pt. Output into neighbours_file which can be viewed from tecplot with OcTreeNeighbours.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: octree.cc:4356
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:342
static int get_the_other_face(const unsigned &n1, const unsigned &n2, const unsigned &nnode1d, const int &face)
If an edge is bordered by the nodes whose local numbers are n1 and n2 in an element with nnode1d node...
Definition: octree.cc:433
OcTree()
Default constructor (empty and broken)
Definition: octree.h:371
static void doc_face_neighbours(Vector< Tree * > forest_nodes_pt, std::ofstream &neighbours_file, std::ofstream &neighbours_txt_file, double &max_error)
Doc/check all face neighbours of octree (nodes) contained in the Vector forest_node_pt. Output into neighbours_file which can be viewed from tecplot with OcTreeNeighbours.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: octree.cc:4099
static DenseMatrix< double > S_directlo
Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son spe...
Definition: octree.h:563
OcTreeRoot(RefineableElement *const &object_pt)
Constructor for the root octree: Pass pointer to the RefineableQElement<3> that is represented by the...
Definition: octree.h:636
static DenseMatrix< double > S_step_edge
Each edge of the RefineableQElement<3> that is represented by the octree is parametrised by one (of the...
Definition: octree.h:583
OcTree(RefineableElement *const &object_pt)
Constructor for empty (root) tree: no father, no sons; just pass a pointer to its object (a Refineabl...
Definition: octree.h:385
OcTree(const OcTree &dummy)
Broken copy constructor.
Definition: octree.h:110
std::map< TreeRoot *, int > Right_equivalent
Map giving the Right equivalent of the neighbour specified by pointer: When viewed from the current o...
Definition: octree.h:629
static void mult_mat_vect(const DenseMatrix< int > &mat, const Vector< int > &vect1, Vector< int > &vect2)
Helper function: Performs the operation : vect2 = mat*vect1.
Definition: octree.cc:533
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:318
OcTree(RefineableElement *const &object_pt, Tree *const &father_pt, const int &son_type)
Constructor for tree that has a father: Pass it the pointer to its object, the pointer to its father ...
Definition: octree.h:394
OcTreeRoot * oc_face_neigh_pt(const unsigned &i, const int &direction)
Given the number i of the root octree in this forest, return pointer to its face neighbour in the spe...
Definition: octree.h:949
std::map< int, TreeRoot * > Neighbour_pt
Map of pointers to the neighbouring TreeRoots: Neighbour_pt[direction] returns the pointer to the Tre...
Definition: tree.h:311
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:241
int son_type() const
Return son type.
Definition: tree.h:209
void set_up_equivalent(TreeRoot *tree_root_pt, const int &dir)
Set up equivalent of the neighbours specified by pointer: When viewed from the current octree's neigh...
Definition: octree.h:791
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
static DenseMatrix< bool > Is_adjacent
Array of direction/octant adjacency scheme: Is_adjacent(direction,octant): Is face/edge direction adj...
Definition: octree.h:519
Vector< TreeRoot * > oc_edge_neigh_pt(const unsigned &i, const int &direction)
Given the number i of the root octree in this forest, return the vector of pointers to the true edge ...
Definition: octree.h:972
OcTreeForest(const OcTreeForest &dummy)
Broken copy constructor.
Definition: octree.h:905
static int node_number_to_vertex(const unsigned &n, const unsigned &nnode1d)
Return the vertex [LDB,RDB,...] of local (vertex) node n.
Definition: octree.cc:375
unsigned self_test()
Self-test: Check all neighbours. Return success (0) if the max. distance between corresponding points...
Definition: octree.cc:4022
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
OcTreeRoot * octree_pt(const unsigned &i) const
Return pointer to i-th OcTree in forest (Performs a dynamic cast from the TreeRoot to a OcTreeRoot)...
Definition: octree.h:942
std::map< TreeRoot *, int > Up_equivalent
Map giving the Up equivalent of the neighbour specified by pointer: When viewed from the current octr...
Definition: octree.h:621
void operator=(const OcTree &)
Broken assignment operator.
Definition: octree.h:116
void add_edge_neighbour_pt(TreeRoot *oc_tree_root_pt, const unsigned &edge_direction)
Add pointer to the edge-neighbouring [Oc]TreeRoot in the (enumerated) (edge) direction – maintains un...
Definition: octree.h:726
OcTreeForest()
Default constructor (empty and broken)
Definition: octree.h:897
int direction_of_neighbour(TreeRoot *octree_root_pt)
If octree_root_pt is a neighbour, return the direction [faces L/R/F/B/U/D or edges DB/UP/...
Definition: octree.h:826
static DenseMatrix< double > S_directhi
Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son spe...
Definition: octree.h:571
static void mult_mat_mat(const DenseMatrix< int > &mat1, const DenseMatrix< int > &mat2, DenseMatrix< int > &mat3)
Helper function: Performs the operation : mat3=mat1*mat2.
Definition: octree.cc:509
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:221
static DenseMatrix< int > Common_face
Determine common face of edges or octants. Slightly bizarre lookup scheme from Samet's book...
Definition: octree.h:527
static DenseMatrix< double > S_direct_edge
Relative to the left/down/back vertex in any (father) octree, the corresponding vertex in the son spe...
Definition: octree.h:590
virtual ~OcTreeForest()
Destructor: Delete the constituent octrees (and thus the associated objects!)
Definition: octree.h:918
static Vector< int > Sini
Entry in rotation matrix sin(i*90)
Definition: octree.h:513
static void construct_rotation_matrix(int &axis, int &angle, DenseMatrix< int > &mat)
This constructs the rotation matrix of the rotation around the axis axis with an angle of angle*90...
Definition: octree.cc:476
Vector< TreeRoot * > edge_neighbour_pt(const unsigned &edge_direction)
Return vector of pointers to the edge-neighbouring TreeRoots in the (enumerated) (edge) direction...
Definition: octree.h:674
int right_equivalent(TreeRoot *tree_root_pt)
The same thing as up_equivalent, but for the right direction: When viewed from the current octree nei...
Definition: octree.h:802
void operator=(const OcTreeForest &)
Broken assignment operator.
Definition: octree.h:911