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