libMesh
tree_node.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 #include <set>
22 
23 // Local includes
24 #include "libmesh/libmesh_config.h"
25 #include "libmesh/tree_node.h"
26 #include "libmesh/mesh_base.h"
27 #include "libmesh/elem.h"
28 
29 namespace libMesh
30 {
31 
32 // ------------------------------------------------------------
33 // TreeNode class methods
34 template <unsigned int N>
35 bool TreeNode<N>::insert (const Node * nd)
36 {
37  libmesh_assert(nd);
38  libmesh_assert_less (nd->id(), mesh.n_nodes());
39 
40  // Return if we don't bound the node
41  if (!this->bounds_node(nd))
42  return false;
43 
44  // Add the node to ourself if we are active
45  if (this->active())
46  {
47  nodes.push_back (nd);
48 
49  // Refine ourself if we reach the target bin size for a TreeNode.
50  if (nodes.size() == tgt_bin_size)
51  this->refine();
52 
53  return true;
54  }
55 
56  // If we are not active simply pass the node along to
57  // our children
58  libmesh_assert_equal_to (children.size(), N);
59 
60  bool was_inserted = false;
61  for (unsigned int c=0; c<N; c++)
62  if (children[c]->insert (nd))
63  was_inserted = true;
64  return was_inserted;
65 }
66 
67 
68 
69 template <unsigned int N>
70 bool TreeNode<N>::insert (const Elem * elem)
71 {
72  libmesh_assert(elem);
73 
74  // We first want to find the corners of the cuboid surrounding the cell.
75  const BoundingBox bbox = elem->loose_bounding_box();
76 
77  // Next, find out whether this cuboid has got non-empty intersection
78  // with the bounding box of the current tree node.
79  //
80  // If not, we should not care about this element.
81  if (!this->bounding_box.intersects(bbox))
82  return false;
83 
84  // Only add the element if we are active
85  if (this->active())
86  {
87  elements.push_back (elem);
88 
89 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
90 
91  // flag indicating this node contains
92  // infinite elements
93  if (elem->infinite())
94  this->contains_ifems = true;
95 
96 #endif
97 
98  unsigned int element_count = elements.size();
100  {
101  const std::set<unsigned char> & elem_dimensions = mesh.elem_dimensions();
102  if (elem_dimensions.size() > 1)
103  {
104  element_count = 0;
105  unsigned char highest_dim_elem = *elem_dimensions.rbegin();
106  for (std::size_t i=0; i<elements.size(); i++)
107  {
108  if (elements[i]->dim() == highest_dim_elem)
109  {
110  element_count++;
111  }
112  }
113  }
114  }
115 
116  // Refine ourself if we reach the target bin size for a TreeNode.
117  if (element_count == tgt_bin_size)
118  this->refine();
119 
120  return true;
121  }
122 
123  // If we are not active simply pass the element along to
124  // our children
125  libmesh_assert_equal_to (children.size(), N);
126 
127  bool was_inserted = false;
128  for (unsigned int c=0; c<N; c++)
129  if (children[c]->insert (elem))
130  was_inserted = true;
131  return was_inserted;
132 }
133 
134 
135 
136 template <unsigned int N>
138 {
139  // Huh? better be active...
140  libmesh_assert (this->active());
141  libmesh_assert (children.empty());
142 
143  // A TreeNode<N> has by definition N children
144  children.resize(N);
145 
146  // Scale up the target bin size in child TreeNodes if we have reached
147  // the maximum number of refinement levels.
148  unsigned int new_target_bin_size = tgt_bin_size;
149  if (level() >= target_bin_size_increase_level)
150  {
151  new_target_bin_size *= 2;
152  }
153 
154  for (unsigned int c=0; c<N; c++)
155  {
156  // Create the child and set its bounding box.
157  children[c] = new TreeNode<N> (mesh, new_target_bin_size, this);
158  children[c]->set_bounding_box(this->create_bounding_box(c));
159 
160  // Pass off our nodes to our children
161  for (std::size_t n=0; n<nodes.size(); n++)
162  children[c]->insert(nodes[n]);
163 
164  // Pass off our elements to our children
165  for (std::size_t e=0; e<elements.size(); e++)
166  children[c]->insert(elements[e]);
167  }
168 
169  // We don't need to store nodes or elements any more, they have been
170  // added to the children. Use the "swap trick" to actually reduce
171  // the capacity of these vectors.
172  std::vector<const Node *>().swap(nodes);
173  std::vector<const Elem *>().swap(elements);
174 
175  libmesh_assert_equal_to (nodes.capacity(), 0);
176  libmesh_assert_equal_to (elements.capacity(), 0);
177 }
178 
179 
180 
181 template <unsigned int N>
182 void TreeNode<N>::set_bounding_box (const std::pair<Point, Point> & bbox)
183 {
184  bounding_box = bbox;
185 }
186 
187 
188 
189 template <unsigned int N>
191  Real relative_tol) const
192 {
193  libmesh_assert(nd);
194  return bounds_point(*nd, relative_tol);
195 }
196 
197 
198 
199 template <unsigned int N>
201  Real relative_tol) const
202 {
203  const Point & min = bounding_box.first;
204  const Point & max = bounding_box.second;
205 
206  const Real tol = (max - min).norm() * relative_tol;
207 
208  if ((p(0) >= min(0) - tol)
209  && (p(0) <= max(0) + tol)
210 #if LIBMESH_DIM > 1
211  && (p(1) >= min(1) - tol)
212  && (p(1) <= max(1) + tol)
213 #endif
214 #if LIBMESH_DIM > 2
215  && (p(2) >= min(2) - tol)
216  && (p(2) <= max(2) + tol)
217 #endif
218  )
219  return true;
220 
221  return false;
222 }
223 
224 
225 
226 template <unsigned int N>
228 TreeNode<N>::create_bounding_box (unsigned int c) const
229 {
230  switch (N)
231  {
232  // How to refine an OctTree Node
233  case 8:
234  {
235  const Real xmin = bounding_box.first(0);
236  const Real ymin = bounding_box.first(1);
237  const Real zmin = bounding_box.first(2);
238 
239  const Real xmax = bounding_box.second(0);
240  const Real ymax = bounding_box.second(1);
241  const Real zmax = bounding_box.second(2);
242 
243  const Real xc = .5*(xmin + xmax);
244  const Real yc = .5*(ymin + ymax);
245  const Real zc = .5*(zmin + zmax);
246 
247  switch (c)
248  {
249  case 0:
250  return BoundingBox (Point(xmin, ymin, zmin),
251  Point(xc, yc, zc));
252  case 1:
253  return BoundingBox (Point(xc, ymin, zmin),
254  Point(xmax, yc, zc));
255  case 2:
256  return BoundingBox (Point(xmin, yc, zmin),
257  Point(xc, ymax, zc));
258  case 3:
259  return BoundingBox (Point(xc, yc, zmin),
260  Point(xmax, ymax, zc));
261  case 4:
262  return BoundingBox (Point(xmin, ymin, zc),
263  Point(xc, yc, zmax));
264  case 5:
265  return BoundingBox (Point(xc, ymin, zc),
266  Point(xmax, yc, zmax));
267  case 6:
268  return BoundingBox (Point(xmin, yc, zc),
269  Point(xc, ymax, zmax));
270  case 7:
271  return BoundingBox (Point(xc, yc, zc),
272  Point(xmax, ymax, zmax));
273  default:
274  libmesh_error_msg("c >= N! : " << c);
275  }
276 
277  break;
278  } // case 8
279 
280  // How to refine an QuadTree Node
281  case 4:
282  {
283  const Real xmin = bounding_box.first(0);
284  const Real ymin = bounding_box.first(1);
285 
286  const Real xmax = bounding_box.second(0);
287  const Real ymax = bounding_box.second(1);
288 
289  const Real xc = .5*(xmin + xmax);
290  const Real yc = .5*(ymin + ymax);
291 
292  switch (c)
293  {
294  case 0:
295  return BoundingBox (Point(xmin, ymin),
296  Point(xc, yc));
297  case 1:
298  return BoundingBox (Point(xc, ymin),
299  Point(xmax, yc));
300  case 2:
301  return BoundingBox (Point(xmin, yc),
302  Point(xc, ymax));
303  case 3:
304  return BoundingBox (Point(xc, yc),
305  Point(xmax, ymax));
306  default:
307  libmesh_error_msg("c >= N!");
308  }
309 
310  break;
311  } // case 4
312 
313  // How to refine a BinaryTree Node
314  case 2:
315  {
316  const Real xmin = bounding_box.first(0);
317 
318  const Real xmax = bounding_box.second(0);
319 
320  const Real xc = .5*(xmin + xmax);
321 
322  switch (c)
323  {
324  case 0:
325  return BoundingBox (Point(xmin),
326  Point(xc));
327  case 1:
328  return BoundingBox (Point(xc),
329  Point(xmax));
330  default:
331  libmesh_error_msg("c >= N!");
332  }
333 
334  break;
335  } // case 2
336 
337  default:
338  libmesh_error_msg("Only implemented for Octrees, QuadTrees, and Binary Trees!");
339  }
340 
341  libmesh_error_msg("We'll never get here!");
342  Point min, max;
343  return BoundingBox (min, max);
344 }
345 
346 
347 
348 template <unsigned int N>
349 void TreeNode<N>::print_nodes(std::ostream & out_stream) const
350 {
351  if (this->active())
352  {
353  out_stream << "TreeNode Level: " << this->level() << std::endl;
354 
355  for (std::size_t n=0; n<nodes.size(); n++)
356  out_stream << " " << nodes[n]->id();
357 
358  out_stream << std::endl << std::endl;
359  }
360  else
361  {
362  for (std::size_t child=0; child<children.size(); child++)
363  children[child]->print_nodes();
364  }
365 }
366 
367 
368 
369 template <unsigned int N>
370 void TreeNode<N>::print_elements(std::ostream & out_stream) const
371 {
372  if (this->active())
373  {
374  out_stream << "TreeNode Level: " << this->level() << std::endl;
375 
376  for (std::vector<const Elem *>::const_iterator pos=elements.begin();
377  pos != elements.end(); ++pos)
378  out_stream << " " << *pos;
379 
380  out_stream << std::endl << std::endl;
381  }
382  else
383  {
384  for (std::size_t child=0; child<children.size(); child++)
385  children[child]->print_elements();
386  }
387 }
388 
389 
390 
391 template <unsigned int N>
392 void TreeNode<N>::transform_nodes_to_elements (std::vector<std::vector<const Elem *>> & nodes_to_elem)
393 {
394  if (this->active())
395  {
396  elements.clear();
397 
398  // Temporarily use a set. Since multiple nodes
399  // will likely map to the same element we use a
400  // set to eliminate the duplication.
401  std::set<const Elem *> elements_set;
402 
403  for (std::size_t n=0; n<nodes.size(); n++)
404  {
405  // the actual global node number we are replacing
406  // with the connected elements
407  const dof_id_type node_number = nodes[n]->id();
408 
409  libmesh_assert_less (node_number, mesh.n_nodes());
410  libmesh_assert_less (node_number, nodes_to_elem.size());
411 
412  for (std::size_t e=0; e<nodes_to_elem[node_number].size(); e++)
413  elements_set.insert(nodes_to_elem[node_number][e]);
414  }
415 
416  // Done with the nodes.
417  std::vector<const Node *>().swap(nodes);
418 
419  // Now the set is built. We can copy this to the
420  // vector. Note that the resulting vector will
421  // already be sorted, and will require less memory
422  // than the set.
423  elements.reserve(elements_set.size());
424 
425  for (std::set<const Elem *>::iterator pos=elements_set.begin();
426  pos != elements_set.end(); ++pos)
427  {
428  elements.push_back(*pos);
429 
430 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
431 
432  // flag indicating this node contains
433  // infinite elements
434  if ((*pos)->infinite())
435  this->contains_ifems = true;
436 
437 #endif
438  }
439  }
440  else
441  {
442  for (std::size_t child=0; child<children.size(); child++)
443  children[child]->transform_nodes_to_elements (nodes_to_elem);
444  }
445 
446 }
447 
448 
449 
450 template <unsigned int N>
451 unsigned int TreeNode<N>::n_active_bins() const
452 {
453  if (this->active())
454  return 1;
455 
456  else
457  {
458  unsigned int sum=0;
459 
460  for (std::size_t c=0; c<children.size(); c++)
461  sum += children[c]->n_active_bins();
462 
463  return sum;
464  }
465 }
466 
467 
468 
469 template <unsigned int N>
470 const Elem *
472  const std::set<subdomain_id_type> * allowed_subdomains,
473  Real relative_tol) const
474 {
475  if (this->active())
476  {
477  // Only check our children if the point is in our bounding box
478  // or if the node contains infinite elements
479  if (this->bounds_point(p, relative_tol) || this->contains_ifems)
480  // Search the active elements in the active TreeNode.
481  for (std::vector<const Elem *>::const_iterator pos=elements.begin();
482  pos != elements.end(); ++pos)
483  if (!allowed_subdomains || allowed_subdomains->count((*pos)->subdomain_id()))
484  if ((*pos)->active() && (*pos)->contains_point(p, relative_tol))
485  return *pos;
486 
487  // The point was not found in any element
488  return libmesh_nullptr;
489  }
490  else
491  return this->find_element_in_children(p,allowed_subdomains,
492  relative_tol);
493 
494  libmesh_error_msg("We'll never get here!");
495  return libmesh_nullptr;
496 }
497 
498 
499 
500 
501 template <unsigned int N>
503  const std::set<subdomain_id_type> * allowed_subdomains,
504  Real relative_tol) const
505 {
506  libmesh_assert (!this->active());
507 
508  searched_child.assign(children.size(), false);
509 
510  // First only look in the children whose bounding box
511  // contain the point p.
512  for (std::size_t c=0; c<children.size(); c++)
513  if (children[c]->bounds_point(p, relative_tol))
514  {
515  const Elem * e =
516  children[c]->find_element(p,allowed_subdomains,
517  relative_tol);
518 
519  if (e != libmesh_nullptr)
520  return e;
521 
522  // If we get here then a child that bounds the
523  // point does not have any elements that contain
524  // the point. So, we will search all our children.
525  // However, we have already searched child c so there
526  // is no use searching her again.
527  searched_child[c] = true;
528  }
529 
530 
531  // If we get here then our child whose bounding box
532  // was searched and did not find any elements containing
533  // the point p. So, let's look at the other children
534  // but exclude the one we have already searched.
535  for (std::size_t c=0; c<children.size(); c++)
536  if (!searched_child[c])
537  {
538  const Elem * e =
539  children[c]->find_element(p,allowed_subdomains,
540  relative_tol);
541 
542  if (e != libmesh_nullptr)
543  return e;
544  }
545 
546  // If we get here we have searched all our children.
547  // Since this process was started at the root node then
548  // we have searched all the elements in the tree without
549  // success. So, we should return NULL since at this point
550  // _no_ elements in the tree claim to contain point p.
551 
552  return libmesh_nullptr;
553 }
554 
555 
556 
557 // ------------------------------------------------------------
558 // Explicit Instantiations
559 template class TreeNode<2>;
560 template class TreeNode<4>;
561 template class TreeNode<8>;
562 
563 } // namespace libMesh
A Node is like a Point, but with more information.
Definition: node.h:52
bool bounds_node(const Node *nd, Real relative_tol=0) const
Definition: tree_node.C:190
void refine()
Refine the tree node into N children if it contains more than tol nodes.
Definition: tree_node.C:137
bool intersects(const BoundingBox &) const
Definition: bounding_box.C:33
void print_nodes(std::ostream &out=libMesh::out) const
Prints the contents of the node_numbers vector if we are active.
Definition: tree_node.C:349
libMesh::BoundingBox create_bounding_box(const MeshBase &mesh)
The same functionality as the deprecated MeshTools::bounding_box().
Definition: mesh_tools.C:329
unsigned int dim
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
long double max(long double a, double b)
libmesh_assert(j)
bool insert(const Node *nd)
Tries to insert Node nd into the TreeNode.
Definition: tree_node.C:35
virtual bool infinite() const =0
BoundingBox bounding_box(const MeshBase &mesh)
Definition: mesh_tools.C:316
unsigned int n_active_bins() const
Definition: tree_node.C:451
void print_elements(std::ostream &out=libMesh::out) const
Prints the contents of the elements set if we are active.
Definition: tree_node.C:370
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:206
bool bounds_point(const Point &p, Real relative_tol=0) const
Definition: tree_node.C:200
This class defines a node on a tree.
Definition: tree_node.h:52
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void transform_nodes_to_elements(std::vector< std::vector< const Elem * >> &nodes_to_elem)
Transforms node numbers to element pointers.
Definition: tree_node.C:392
bool get_count_lower_dim_elems_in_point_locator() const
Get the current value of _count_lower_dim_elems_in_point_locator.
Definition: mesh_base.C:569
X_input swap(X_system)
const Elem * find_element_in_children(const Point &p, const std::set< subdomain_id_type > *allowed_subdomains, Real relative_tol) const
Look for point p in our children, optionally restricted to a set of allowed subdomains.
Definition: tree_node.C:502
void set_bounding_box(const std::pair< Point, Point > &bbox)
Sets the bounding box;.
Definition: tree_node.C:182
dof_id_type id() const
Definition: dof_object.h:632
virtual dof_id_type n_nodes() const =0
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
BoundingBox create_bounding_box(unsigned int c) const
Constructs the bounding box for child c.
Definition: tree_node.C:228
const Elem * find_element(const Point &p, const std::set< subdomain_id_type > *allowed_subdomains=libmesh_nullptr, Real relative_tol=TOLERANCE) const
Definition: tree_node.C:471
virtual BoundingBox loose_bounding_box() const
Definition: elem.C:2997
uint8_t dof_id_type
Definition: id_types.h:64