tree_node.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 void 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;
43 
44  // Only add the node 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 
54  // If we are not active simply pass the node along to
55  // our children
56  else
57  {
58  libmesh_assert_equal_to (children.size(), N);
59 
60  for (unsigned int c=0; c<N; c++)
61  children[c]->insert (nd);
62  }
63 }
64 
65 
66 
67 template <unsigned int N>
68 void TreeNode<N>::insert (const Elem* elem)
69 {
70  libmesh_assert(elem);
71 
72  /* We first want to find the corners of the cuboid surrounding the
73  cell. */
74  Point minCoord = elem->point(0);
75  Point maxCoord = minCoord;
76  unsigned int dim = elem->dim();
77  for(unsigned int i=elem->n_nodes()-1; i>0; i--)
78  {
79  Point p = elem->point(i);
80  for(unsigned int d=0; d<dim; d++)
81  {
82  if(minCoord(d)>p(d)) minCoord(d) = p(d);
83  if(maxCoord(d)<p(d)) maxCoord(d) = p(d);
84  }
85  }
86 
87  /* Next, find out whether this cuboid has got non-empty intersection
88  with the bounding box of the current tree node. */
89  bool intersects = true;
90  for(unsigned int d=0; d<dim; d++)
91  {
92  if(maxCoord(d)<this->bounding_box.first(d) ||
93  minCoord(d)>this->bounding_box.second(d))
94  {
95  intersects = false;
96  }
97  }
98 
99  /* If not, we should not care about this element. */
100  if(!intersects)
101  {
102  return;
103  }
104 
105  // Only add the element if we are active
106  if (this->active())
107  {
108  elements.push_back (elem);
109 
110 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
111 
112  // flag indicating this node contains
113  // infinite elements
114  if (elem->infinite())
115  this->contains_ifems = true;
116 
117 #endif
118 
119  // Refine ourself if we reach the target bin size for a TreeNode.
120  if (elements.size() == tgt_bin_size)
121  this->refine();
122  }
123 
124  // If we are not active simply pass the element along to
125  // our children
126  else
127  {
128  libmesh_assert_equal_to (children.size(), N);
129 
130  for (unsigned int c=0; c<N; c++)
131  children[c]->insert (elem);
132  }
133 }
134 
135 
136 
137 template <unsigned int N>
139 {
140  // Huh? better be active...
141  libmesh_assert (this->active());
142  libmesh_assert (children.empty());
143 
144  // A TreeNode<N> has by definition N children
145  children.resize(N);
146 
147  for (unsigned int c=0; c<N; c++)
148  {
149  // Create the child and set its bounding box.
150  children[c] = new TreeNode<N> (mesh, tgt_bin_size, this);
151  children[c]->set_bounding_box(this->create_bounding_box(c));
152 
153  // Pass off our nodes to our children
154  for (unsigned int n=0; n<nodes.size(); n++)
155  children[c]->insert(nodes[n]);
156 
157  // Pass off our elements to our children
158  for (unsigned int e=0; e<elements.size(); e++)
159  children[c]->insert(elements[e]);
160  }
161 
162  // We don't need to store nodes or elements any more,
163  // they have been added to the children.
164  // Note that we cannot use std::vector<>::clear() here
165  // since that in general does not reduce capacity!!
166  // That would be a bad, bad thing.
167  std::vector<const Node*>().swap(nodes);
168  std::vector<const Elem*>().swap(elements);
169 
170  libmesh_assert_equal_to (nodes.capacity(), 0);
171  libmesh_assert_equal_to (elements.capacity(), 0);
172 }
173 
174 
175 
176 template <unsigned int N>
177 void TreeNode<N>::set_bounding_box (const std::pair<Point, Point>& bbox)
178 {
179  bounding_box = bbox;
180 }
181 
182 
183 
184 template <unsigned int N>
185 bool TreeNode<N>::bounds_point (const Point& p) const
186 {
187  const Point& min = bounding_box.first;
188  const Point& max = bounding_box.second;
189 
190 
191  if ((p(0) >= min(0))
192  && (p(0) <= max(0))
193 #if LIBMESH_DIM > 1
194  && (p(1) >= min(1))
195  && (p(1) <= max(1))
196 #endif
197 #if LIBMESH_DIM > 2
198  && (p(2) >= min(2))
199  && (p(2) <= max(2))
200 #endif
201  )
202  return true;
203 
204  return false;
205 }
206 
207 
208 
209 template <unsigned int N>
210 std::pair<Point, Point>
211 TreeNode<N>::create_bounding_box (const unsigned int c) const
212 {
213  switch (N)
214  {
215 
216  // How to refine an OctTree Node
217  case 8:
218  {
219  const Real xmin = bounding_box.first(0);
220  const Real ymin = bounding_box.first(1);
221  const Real zmin = bounding_box.first(2);
222 
223  const Real xmax = bounding_box.second(0);
224  const Real ymax = bounding_box.second(1);
225  const Real zmax = bounding_box.second(2);
226 
227  const Real xc = xmin + .5*(xmax - xmin);
228  const Real yc = ymin + .5*(ymax - ymin);
229  const Real zc = zmin + .5*(zmax - zmin);
230 
231 
232  switch (c)
233  {
234 
235  case 0:
236  {
237  const Point min(xmin, ymin, zmin);
238  const Point max(xc, yc, zc);
239  return std::make_pair (min, max);
240  break;
241  }
242 
243  case 1:
244  {
245  const Point min(xc, ymin, zmin);
246  const Point max(xmax, yc, zc);
247  return std::make_pair (min, max);
248  break;
249  }
250 
251  case 2:
252  {
253  const Point min(xmin, yc, zmin);
254  const Point max(xc, ymax, zc);
255  return std::make_pair (min, max);
256  break;
257  }
258 
259  case 3:
260  {
261  const Point min(xc, yc, zmin);
262  const Point max(xmax, ymax, zc);
263  return std::make_pair (min, max);
264  break;
265  }
266 
267  case 4:
268  {
269  const Point min(xmin, ymin, zc);
270  const Point max(xc, yc, zmax);
271  return std::make_pair (min, max);
272  break;
273  }
274 
275  case 5:
276  {
277  const Point min(xc, ymin, zc);
278  const Point max(xmax, yc, zmax);
279  return std::make_pair (min, max);
280  break;
281  }
282 
283  case 6:
284  {
285  const Point min(xmin, yc, zc);
286  const Point max(xc, ymax, zmax);
287  return std::make_pair (min, max);
288  break;
289  }
290 
291  case 7:
292  {
293  const Point min(xc, yc, zc);
294  const Point max(xmax, ymax, zmax);
295  return std::make_pair (min, max);
296  break;
297  }
298 
299  default:
300  libMesh::err << "c >= N! : " << c
301  << std::endl;
302  libmesh_error();
303  }
304 
305 
306 
307  break;
308  } // case 8
309 
310  // How to refine an QuadTree Node
311  case 4:
312  {
313  const Real xmin = bounding_box.first(0);
314  const Real ymin = bounding_box.first(1);
315 
316  const Real xmax = bounding_box.second(0);
317  const Real ymax = bounding_box.second(1);
318 
319  const Real xc = xmin + .5*(xmax - xmin);
320  const Real yc = ymin + .5*(ymax - ymin);
321 
322  switch (c)
323  {
324  case 0:
325  {
326  const Point min(xmin, ymin);
327  const Point max(xc, yc);
328  return std::make_pair (min, max);
329  break;
330  }
331 
332  case 1:
333  {
334  const Point min(xc, ymin);
335  const Point max(xmax, yc);
336  return std::make_pair (min, max);
337  break;
338  }
339 
340  case 2:
341  {
342  const Point min(xmin, yc);
343  const Point max(xc, ymax);
344  return std::make_pair (min, max);
345  break;
346  }
347 
348  case 3:
349  {
350  const Point min(xc, yc);
351  const Point max(xmax, ymax);
352  return std::make_pair (min, max);
353  break;
354  }
355 
356  default:
357  libMesh::err << "c >= N!" << std::endl;
358  libmesh_error();
359 
360  }
361 
362  break;
363  } // case 4
364 
365  // How to refine a BinaryTree Node
366  case 2:
367  {
368  const Real xmin = bounding_box.first(0);
369 
370  const Real xmax = bounding_box.second(0);
371 
372  const Real xc = xmin + .5*(xmax - xmin);
373 
374  switch (c)
375  {
376  case 0:
377  {
378  return std::make_pair (Point(xmin), Point(xc));
379  break;
380  }
381 
382  case 1:
383  {
384  return std::make_pair (Point(xc), Point(xmax));
385  break;
386  }
387 
388  default:
389  libMesh::err << "c >= N!" << std::endl;
390  libmesh_error();
391  }
392 
393  break;
394  } // case 2
395 
396 
397  default:
398  libMesh::err << "Only implemented for Octrees, QuadTrees, and Binary Trees!" << std::endl;
399  libmesh_error();
400 
401  }
402 
403  // How did we get here?
404  libmesh_error();
405 
406  Point min, max;
407  return std::make_pair (min, max);
408 }
409 
410 
411 
412 template <unsigned int N>
413 void TreeNode<N>::print_nodes(std::ostream& out_stream) const
414 {
415  if (this->active())
416  {
417  out_stream << "TreeNode Level: " << this->level() << std::endl;
418 
419  for (unsigned int n=0; n<nodes.size(); n++)
420  out_stream << " " << nodes[n]->id();
421 
422  out_stream << std::endl << std::endl;
423 
424  }
425  else
426  {
427  for (unsigned int child=0; child<children.size(); child++)
428  children[child]->print_nodes();
429  }
430 }
431 
432 
433 
434 template <unsigned int N>
435 void TreeNode<N>::print_elements(std::ostream& out_stream) const
436 {
437  if (this->active())
438  {
439  out_stream << "TreeNode Level: " << this->level() << std::endl;
440 
441  for (std::vector<const Elem*>::const_iterator pos=elements.begin();
442  pos != elements.end(); ++pos)
443  out_stream << " " << *pos;
444 
445  out_stream << std::endl << std::endl;
446  }
447  else
448  {
449  for (unsigned int child=0; child<children.size(); child++)
450  children[child]->print_elements();
451  }
452 }
453 
454 
455 
456 template <unsigned int N>
457 void TreeNode<N>::transform_nodes_to_elements (std::vector<std::vector<const Elem*> >& nodes_to_elem)
458 {
459  if (this->active())
460  {
461  elements.clear();
462 
463  // Temporarily use a set. Since multiple nodes
464  // will likely map to the same element we use a
465  // set to eliminate the duplication.
466  std::set<const Elem*> elements_set;
467 
468  for (unsigned int n=0; n<nodes.size(); n++)
469  {
470  // the actual global node number we are replacing
471  // with the connected elements
472  const dof_id_type node_number = nodes[n]->id();
473 
474  libmesh_assert_less (node_number, mesh.n_nodes());
475  libmesh_assert_less (node_number, nodes_to_elem.size());
476 
477  for (unsigned int e=0; e<nodes_to_elem[node_number].size(); e++)
478  elements_set.insert(nodes_to_elem[node_number][e]);
479  }
480 
481  // Done with the nodes.
482  std::vector<const Node*>().swap(nodes);
483 
484  // Now the set is built. We can copy this to the
485  // vector. Note that the resulting vector will
486  // already be sorted, and will require less memory
487  // than the set.
488  elements.reserve(elements_set.size());
489 
490  for (std::set<const Elem*>::iterator pos=elements_set.begin();
491  pos != elements_set.end(); ++pos)
492  {
493  elements.push_back(*pos);
494 
495 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
496 
497  // flag indicating this node contains
498  // infinite elements
499  if ((*pos)->infinite())
500  this->contains_ifems = true;
501 
502 #endif
503  }
504  }
505  else
506  {
507  for (unsigned int child=0; child<children.size(); child++)
508  children[child]->transform_nodes_to_elements (nodes_to_elem);
509  }
510 
511 }
512 
513 
514 
515 template <unsigned int N>
516 unsigned int TreeNode<N>::n_active_bins() const
517 {
518  if (this->active())
519  return 1;
520 
521  else
522  {
523  unsigned int sum=0;
524 
525  for (unsigned int c=0; c<children.size(); c++)
526  sum += children[c]->n_active_bins();
527 
528  return sum;
529  }
530 }
531 
532 
533 
534 template <unsigned int N>
535 const Elem* TreeNode<N>::find_element(const Point& p) const
536 {
537  if (this->active())
538  {
539  // Only check our children if the point is in our bounding box
540  // or if the node contains infinite elements
541  if (this->bounds_point(p) || this->contains_ifems)
542  // Search the active elements in the active TreeNode.
543  for (std::vector<const Elem*>::const_iterator pos=elements.begin();
544  pos != elements.end(); ++pos)
545  if ((*pos)->active())
546  if ((*pos)->contains_point(p))
547  return *pos;
548 
549  // The point was not found in any element
550  return NULL;
551  }
552  else
553  return this->find_element_in_children(p);
554 
555 
556 
557  // Should never get here. See if-else structure
558  // above with return statements that must get executed.
559  libmesh_error();
560 
561  return NULL;
562 }
563 
564 
565 
566 
567 template <unsigned int N>
569 {
570  libmesh_assert (!this->active());
571 
572  unsigned int excluded_child = libMesh::invalid_uint;
573 
574  // First only look in the children whose bounding box
575  // contain the point p. Note that only one child will
576  // bound the point since the bounding boxes are not
577  // overlapping
578  for (unsigned int c=0; c<children.size(); c++)
579  if (children[c]->bounds_point(p))
580  {
581  if (children[c]->active())
582  {
583  const Elem* e = children[c]->find_element(p);
584 
585  if (e != NULL)
586  return e;
587  }
588  else
589  {
590  const Elem* e = children[c]->find_element_in_children(p);
591 
592  if (e != NULL)
593  return e;
594  }
595 
596  // If we get here than the child that bounds the
597  // point does not have any elements that contain
598  // the point. So, we will search all our children.
599  // However, we have already searched child c so there
600  // is no use searching her again.
601  excluded_child = c;
602  }
603 
604 
605  // If we get here then our child whose bounding box
606  // was searched and did not find any elements containing
607  // the point p. So, let's look at the other children
608  // but exclude the one we have already searched.
609  for (unsigned int c=0; c<children.size(); c++)
610  if (c != excluded_child)
611  {
612  if (children[c]->active())
613  {
614  const Elem* e = children[c]->find_element(p);
615 
616  if (e != NULL)
617  return e;
618  }
619  else
620  {
621  const Elem* e = children[c]->find_element_in_children(p);
622 
623  if (e != NULL)
624  return e;
625  }
626  }
627 
628  // If we get here we have searched all our children.
629  // Since this process was started at the root node then
630  // we have searched all the elements in the tree without
631  // success. So, we should return NULL since at this point
632  // _no_ elements in the tree claim to contain point p.
633 
634  return NULL;
635 }
636 
637 
638 
639 // ------------------------------------------------------------
640 // Explicit Instantiations
641 template class TreeNode<2>;
642 template class TreeNode<4>;
643 template class TreeNode<8>;
644 
645 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo