Neko  0.9.0
A portable framework for high-order spectral element flow simulations
aabb_tree.f90
Go to the documentation of this file.
1 ! Copyright (c) 2024, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
33 ! ============================================================================ !
34 ! Original C++ Implementation from:
35 ! https://github.com/JamesRandall/SimpleVoxelEngine/blob/master/voxelEngine/src/AABBTree.h
36 !
37 ! Translated to Fortran by:
38 ! @author Tim Felle Olsen
39 ! @date 9 Feb 2024
40 !
41 ! C++ Code License:
42 ! The MIT License (MIT)
43 !
44 ! Copyright (c) 2017 James Randall
45 !
46 ! Permission is hereby granted, free of charge, to any person obtaining a copy of
47 ! this software and associated documentation files (the "Software"), to deal in
48 ! the Software without restriction, including without limitation the rights to
49 ! use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
50 ! the Software, and to permit persons to whom the Software is furnished to do so,
51 ! subject to the following conditions:
52 !
53 ! The above copyright notice and this permission notice shall be included in all
54 ! copies or substantial portions of the Software.
55 !
56 ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
57 ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
58 ! FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
59 ! COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
60 ! IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
61 ! CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
62 ! ============================================================================ !
63 
70 module aabb_tree
71  use aabb
72  use tri, only: tri_t
73  use num_types, only: rp, dp
74 
75  implicit none
76  private
77 
78  integer, parameter, public :: aabb_null_node = -1
79 
80  ! ========================================================================== !
81  ! Type definitions
82  ! ========================================================================== !
83 
85  type, public :: aabb_node_t
86  private
87  type(aabb_t), public :: aabb
88  integer :: object_index = -1
89 
90  ! tree links
91  integer :: parent_node_index = aabb_null_node
92  integer :: left_node_index = aabb_null_node
93  integer :: right_node_index = aabb_null_node
94 
95  ! node linked list link
96  integer :: next_node_index = aabb_null_node
97 
98  contains
99  procedure, pass(this), public :: init => aabb_node_init
100 
101  ! Getters
102  procedure, pass(this), public :: get_aabb => aabb_node_get_aabb
103  procedure, pass(this), public :: get_object_index => &
105  procedure, pass(this), public :: get_parent_index => &
107  procedure, pass(this), public :: get_left_index => &
109  procedure, pass(this), public :: get_right_index => &
111 
112  ! Unary operations
113  procedure, pass(this), public :: min_distance => aabb_node_min_distance
114 
115  ! Boolean operators
116  procedure, pass(this), public :: is_leaf => aabb_node_is_leaf
117  procedure, pass(this), public :: is_valid => aabb_node_is_valid
118 
119  ! Comparison operators
120  generic :: operator(.lt.) => less
121  generic :: operator(.gt.) => greater
122 
123  procedure, pass(this) :: less => aabb_node_less
124  procedure, pass(this) :: greater => aabb_node_greater
125 
126  end type aabb_node_t
127 
129  type, public :: aabb_tree_t
130  private
131  type(aabb_node_t), allocatable :: nodes(:)
132  integer :: root_node_index = aabb_null_node
133  integer :: allocated_node_count = 0
134  integer :: next_free_node_index = aabb_null_node
135  integer :: node_capacity = 0
136  integer :: growth_size = 1
137 
138  contains
139 
140  ! Initializers
141  procedure, pass(this), public :: init => aabb_tree_init
142  procedure, pass(this), public :: build => aabb_tree_build_tree
143  procedure, pass(this), public :: insert_object => &
145 
146  ! Getters
147  procedure, pass(this), public :: get_size => aabb_tree_get_size
148 
149  procedure, pass(this), public :: get_root_index => &
151  procedure, pass(this), public :: get_parent_index => &
153  procedure, pass(this), public :: get_left_index => &
155  procedure, pass(this), public :: get_right_index => &
157 
158  procedure, pass(this), public :: get_node => aabb_tree_get_node
159  procedure, pass(this), public :: get_root_node => &
161  procedure, pass(this), public :: get_parent_node => &
163  procedure, pass(this), public :: get_left_node => &
165  procedure, pass(this), public :: get_right_node => &
167 
168  procedure, pass(this), public :: get_aabb => aabb_tree_get_aabb
169 
170  procedure, pass(this), public :: query_overlaps => &
172 
173  procedure, pass(this), public :: print => aabb_tree_print
174 
175  ! ----------------------------------------------------------------------- !
176  ! Internal methods
177 
178  procedure, pass(this) :: allocate_node => aabb_tree_allocate_node
179  procedure, pass(this) :: deallocate_node => aabb_tree_deallocate_node
180  procedure, pass(this) :: resize_node_pool => aabb_tree_resize_node_pool
181  procedure, pass(this) :: insert_leaf => aabb_tree_insert_leaf
182 
183  procedure, pass(this) :: fix_upwards_tree => aabb_tree_fix_upwards_tree
184 
185  procedure, pass(this) :: valid_tree => aabb_tree_valid_tree
186 
187  end type aabb_tree_t
188 
189 contains
190 
191  ! ========================================================================== !
192  ! Definitions of node methods
193  ! ========================================================================== !
194 
196  subroutine aabb_node_init(this)
197  class(aabb_node_t), intent(inout) :: this
198 
199  this%object_index = -1
200  this%parent_node_index = aabb_null_node
201  this%left_node_index = aabb_null_node
202  this%right_node_index = aabb_null_node
203  this%next_node_index = aabb_null_node
204  end subroutine aabb_node_init
205 
206  ! -------------------------------------------------------------------------- !
207  ! Getters
208 
210  pure function aabb_node_get_aabb(this) result(res)
211  class(aabb_node_t), intent(in) :: this
212  type(aabb_t) :: res
213 
214  res = this%aabb
215  end function aabb_node_get_aabb
216 
218  pure function aabb_node_get_object_index(this) result(object_index)
219  class(aabb_node_t), intent(in) :: this
220  integer :: object_index
221 
222  object_index = this%object_index
223  end function aabb_node_get_object_index
224 
226  pure function aabb_node_get_parent_index(this) result(parent_index)
227  class(aabb_node_t), intent(in) :: this
228  integer :: parent_index
229 
230  parent_index = this%parent_node_index
231  end function aabb_node_get_parent_index
232 
234  pure function aabb_node_get_left_index(this) result(left_index)
235  class(aabb_node_t), intent(in) :: this
236  integer :: left_index
237 
238  left_index = this%left_node_index
239  end function aabb_node_get_left_index
240 
242  pure function aabb_node_get_right_index(this) result(right_index)
243  class(aabb_node_t), intent(in) :: this
244  integer :: right_index
245 
246  right_index = this%right_node_index
247  end function aabb_node_get_right_index
248 
250  function aabb_node_min_distance(this, p) result(distance)
251  class(aabb_node_t), intent(in) :: this
252  real(kind=dp), dimension(3), intent(in) :: p
253  real(kind=dp) :: distance
254 
255  distance = 0.5_rp * this%aabb%get_diameter() &
256  - norm2(this%aabb%get_center() - p)
257  end function aabb_node_min_distance
258 
259  ! -------------------------------------------------------------------------- !
260  ! Boolean operators
261 
263  pure function aabb_node_is_leaf(this) result(res)
264  class(aabb_node_t), intent(in) :: this
265  logical :: res
266 
267  res = this%left_node_index == aabb_null_node .and. &
268  this%right_node_index == aabb_null_node
269  end function aabb_node_is_leaf
270 
272  pure function aabb_node_is_valid(this) result(valid)
273  class(aabb_node_t), intent(in) :: this
274  logical :: valid
275 
276  if (this%is_leaf()) then
277  valid = &
278  & this%left_node_index .eq. aabb_null_node .and. &
279  & this%right_node_index .eq. aabb_null_node .and. &
280  & this%object_index .gt. 0
281  else
282  valid = &
283  & this%left_node_index .ne. aabb_null_node .and. &
284  & this%right_node_index .ne. aabb_null_node .and. &
285  & this%object_index .eq. -1
286  end if
287 
288  end function aabb_node_is_valid
289 
290  ! -------------------------------------------------------------------------- !
291  ! Comparison operators
292 
294  pure function aabb_node_less(this, other) result(res)
295  class(aabb_node_t), intent(in) :: this
296  class(aabb_node_t), intent(in) :: other
297  logical :: res
298 
299  res = this%aabb .lt. other%aabb
300 
301  end function aabb_node_less
302 
304  pure function aabb_node_greater(this, other) result(res)
305  class(aabb_node_t), intent(in) :: this
306  class(aabb_node_t), intent(in) :: other
307  logical :: res
308 
309  res = this%aabb .gt. other%aabb
310 
311  end function aabb_node_greater
312 
313  ! ========================================================================== !
314  ! Definitions of tree methods
315  ! ========================================================================== !
316 
318  subroutine aabb_tree_init(this, initial_capacity)
319  class(aabb_tree_t), intent(inout) :: this
320  integer, intent(in) :: initial_capacity
321 
322  integer :: i
323 
324  this%root_node_index = aabb_null_node
325  this%allocated_node_count = 0
326  this%next_free_node_index = 1
327  this%node_capacity = initial_capacity
328  this%growth_size = initial_capacity
329 
330  if (allocated(this%nodes)) deallocate(this%nodes)
331  allocate(this%nodes(initial_capacity))
332 
333  do i = 1, initial_capacity
334  this%nodes(i)%next_node_index = i + 1
335  end do
336  this%nodes(initial_capacity)%next_node_index = aabb_null_node
337  end subroutine aabb_tree_init
338 
340  subroutine aabb_tree_build_tree(this, objects)
341  use utils, only: neko_error
342  implicit none
343 
344  class(aabb_tree_t), intent(inout) :: this
345  class(*), dimension(:), intent(in) :: objects
346 
347  integer :: i_obj, i_node, i
348  logical :: done
349 
350  integer :: start_layer, end_layer
351 
352  type(aabb_t), dimension(:), allocatable :: box_list
353  integer, dimension(:), allocatable :: sorted_indices
354 
355  call this%init(size(objects) * 2)
356 
357  ! ------------------------------------------------------------------------ !
358  ! Start by sorting the list of objects, then build a balanced binary tree
359  ! from the sorted list
360 
361  allocate(box_list(size(objects)))
362 
363  do i_obj = 1, size(objects)
364  box_list(i_obj) = get_aabb(objects(i_obj))
365  end do
366  sorted_indices = sort(box_list)
367 
368  do i = 1, size(sorted_indices)
369  i_obj = sorted_indices(i)
370  i_node = this%allocate_node()
371  this%nodes(i_node)%aabb = get_aabb(objects(i_obj))
372  this%nodes(i_node)%object_index = i_obj
373  end do
374 
375 
376  start_layer = 1
377  end_layer = size(objects)
378  done = .false.
379  do while (.not. done)
380 
381  ! build the next layer
382  do i = start_layer, end_layer - 1, 2
383  i_node = this%allocate_node()
384 
385  this%nodes(i_node)%aabb = merge(this%nodes(i)%aabb, &
386  this%nodes(i + 1)%aabb)
387 
388  this%nodes(i_node)%left_node_index = i
389  this%nodes(i_node)%right_node_index = i + 1
390 
391  this%nodes(i)%parent_node_index = i_node
392  this%nodes(i + 1)%parent_node_index = i_node
393  end do
394 
395  ! if the number of nodes is odd, we need to create a new node to hold the
396  ! last node
397  if (mod(end_layer - start_layer, 2) .eq. 0) then
398  i_node = this%allocate_node()
399  this%nodes(i_node)%aabb = this%nodes(end_layer)%aabb
400  this%nodes(i_node)%left_node_index = end_layer
401  this%nodes(i_node)%right_node_index = aabb_null_node
402 
403  this%nodes(end_layer)%parent_node_index = i_node
404  end if
405 
406  ! move to the next layer
407  start_layer = end_layer + 1
408  end_layer = this%allocated_node_count
409 
410  ! If there is only one node left, we are done
411  done = start_layer .eq. end_layer
412  end do
413 
414  ! The last node allocated is the root node
415  this%root_node_index = this%allocated_node_count
416 
417  if (this%get_size() .ne. size(objects)) then
418  print *, "this%get_size() = ", this%get_size()
419  print *, "size(objects) = ", size(objects)
420  call neko_error("Invalid tree size")
421  end if
422 
423  end subroutine aabb_tree_build_tree
424 
425  function sort(array) result(indices)
426  type(aabb_t), dimension(:), intent(in) :: array
427  integer, dimension(:), allocatable :: indices
428  logical, dimension(:), allocatable :: visited
429 
430  integer :: i, imin
431  integer :: minidx
432 
433  allocate(indices(size(array)))
434  allocate(visited(size(array)))
435 
436  visited = .false.
437  indices = 0
438  do i = 1, size(array)
439  minidx = -1
440  do imin = 1, size(array)
441  if (.not. visited(imin) .and. minidx .eq. -1) minidx = imin
442 
443  if (visited(imin) .and. array(imin) .lt. array(minidx)) minidx = imin
444  end do
445 
446  indices(i) = minidx
447  visited(minidx) = .true.
448  end do
449 
450  end function sort
451 
452  ! -------------------------------------------------------------------------- !
453  ! Getters
454 
456  function aabb_tree_get_size(this) result(size)
457  use stack, only: stack_i4_t
458  use utils, only: neko_error
459  class(aabb_tree_t), intent(in) :: this
460  integer :: size
461 
462  type(stack_i4_t) :: simple_stack
463  integer :: idx, tmp
464 
465  call simple_stack%init(this%allocated_node_count)
466  size = 0
467  tmp = this%get_root_index()
468  if (tmp .ne. aabb_null_node) then
469  call simple_stack%push(tmp)
470  end if
471 
472  do while (.not. simple_stack%is_empty())
473  idx = simple_stack%pop()
474  if (idx .eq. aabb_null_node) cycle
475 
476  if (this%nodes(idx)%is_leaf()) then
477  size = size + 1
478  else
479  tmp = this%get_left_index(idx)
480  call simple_stack%push(tmp)
481  tmp = this%get_right_index(idx)
482  call simple_stack%push(tmp)
483  end if
484  end do
485 
486  end function aabb_tree_get_size
487 
488  ! -------------------------------------------------------------------------- !
489  ! Get index of nodes
490 
492  pure function aabb_tree_get_root_index(this) result(root_index)
493  class(aabb_tree_t), intent(in) :: this
494  integer :: root_index
495 
496  root_index = this%root_node_index
497  end function aabb_tree_get_root_index
498 
500  pure function aabb_tree_get_parent_index(this, node_index) &
501  result(parent_index)
502  class(aabb_tree_t), intent(in) :: this
503  integer, intent(in) :: node_index
504  integer :: parent_index
505 
506  parent_index = this%nodes(node_index)%parent_node_index
507  end function aabb_tree_get_parent_index
508 
510  pure function aabb_tree_get_left_index(this, node_index) &
511  result(left_index)
512  class(aabb_tree_t), intent(in) :: this
513  integer, intent(in) :: node_index
514  integer :: left_index
515 
516  left_index = this%nodes(node_index)%left_node_index
517  end function aabb_tree_get_left_index
518 
520  pure function aabb_tree_get_right_index(this, node_index) &
521  result(right_index)
522  class(aabb_tree_t), intent(in) :: this
523  integer, intent(in) :: node_index
524  integer :: right_index
525 
526  right_index = this%nodes(node_index)%right_node_index
527  end function aabb_tree_get_right_index
528 
529  ! -------------------------------------------------------------------------- !
530  ! Get nodes
531 
533  pure function aabb_tree_get_node(this, node_index) result(node)
534  class(aabb_tree_t), intent(in) :: this
535  integer, intent(in) :: node_index
536  type(aabb_node_t) :: node
537 
538  node = this%nodes(node_index)
539  end function aabb_tree_get_node
540 
542  pure function aabb_tree_get_root_node(this) result(root_node)
543  class(aabb_tree_t), intent(in) :: this
544  type(aabb_node_t) :: root_node
545 
546  root_node = this%nodes(this%root_node_index)
547  end function aabb_tree_get_root_node
548 
550  pure function aabb_tree_get_parent_node(this, node_index) &
551  result(parent_node)
552  class(aabb_tree_t), intent(in) :: this
553  integer, intent(in) :: node_index
554  type(aabb_node_t) :: parent_node
555 
556  parent_node = this%nodes(this%nodes(node_index)%parent_node_index)
557  end function aabb_tree_get_parent_node
558 
560  pure function aabb_tree_get_left_node(this, node_index) result(left_node)
561  class(aabb_tree_t), intent(in) :: this
562  integer, intent(in) :: node_index
563  type(aabb_node_t) :: left_node
564 
565  left_node = this%nodes(this%nodes(node_index)%left_node_index)
566  end function aabb_tree_get_left_node
567 
569  pure function aabb_tree_get_right_node(this, node_index) &
570  result(right_node)
571  class(aabb_tree_t), intent(in) :: this
572  integer, intent(in) :: node_index
573  type(aabb_node_t) :: right_node
574 
575  right_node = this%nodes(this%nodes(node_index)%right_node_index)
576  end function aabb_tree_get_right_node
577 
578  pure function aabb_tree_get_aabb(this, node_index) result(out_box)
579  class(aabb_tree_t), intent(in) :: this
580  integer, intent(in) :: node_index
581  type(aabb_t) :: out_box
582 
583  out_box = this%nodes(node_index)%aabb
584  end function aabb_tree_get_aabb
585 
586  ! -------------------------------------------------------------------------- !
587 
589  subroutine aabb_tree_insert_object(this, object, object_index)
590  class(aabb_tree_t), intent(inout) :: this
591  class(*), intent(in) :: object
592  integer, intent(in) :: object_index
593 
594  integer :: node_index
595 
596  node_index = this%allocate_node()
597  this%nodes(node_index)%aabb = get_aabb(object)
598  this%nodes(node_index)%object_index = object_index
599 
600  call this%insert_leaf(node_index)
601  end subroutine aabb_tree_insert_object
602 
604  subroutine aabb_tree_query_overlaps(this, object, object_index, overlaps)
605  use stack, only: stack_i4_t
606  implicit none
607 
608  class(aabb_tree_t), intent(in) :: this
609  class(*), intent(in) :: object
610  integer, intent(in) :: object_index
611  integer, intent(out) :: overlaps(:)
612 
613  type(stack_i4_t) :: simple_stack
614  type(aabb_t) :: object_box
615 
616  integer :: root_index, left_index, right_index
617 
618  integer :: node_index
619 
620  object_box = get_aabb(object)
621  root_index = this%get_root_index()
622 
623  call simple_stack%init()
624  call simple_stack%push(root_index)
625 
626  do while (.not. simple_stack%is_empty())
627  node_index = simple_stack%pop()
628 
629  if (node_index == aabb_null_node) cycle
630 
631  if (this%nodes(node_index)%aabb%overlaps(object_box)) then
632  if (this%nodes(node_index)%is_leaf()) then
633  if (this%nodes(node_index)%object_index .ne. object_index) then
634  overlaps = [this%nodes(node_index)%object_index, overlaps]
635  end if
636  else
637  left_index = this%get_left_index(node_index)
638  call simple_stack%push(left_index)
639  right_index = this%get_right_index(node_index)
640  call simple_stack%push(right_index)
641  end if
642  end if
643  end do
644  end subroutine aabb_tree_query_overlaps
645 
646  ! -------------------------------------------------------------------------- !
647  ! Internal methods
648 
650  function aabb_tree_allocate_node(this) result(node_index)
651  class(aabb_tree_t), intent(inout) :: this
652  integer :: node_index
653 
654  if (this%next_free_node_index == aabb_null_node) then
655  call this%resize_node_pool(this%node_capacity + this%growth_size)
656  end if
657 
658  node_index = this%next_free_node_index
659 
660  associate(new_node => this%nodes(node_index))
661  this%next_free_node_index = new_node%next_node_index
662 
663  new_node%parent_node_index = aabb_null_node
664  new_node%left_node_index = aabb_null_node
665  new_node%right_node_index = aabb_null_node
666 
667  this%next_free_node_index = new_node%next_node_index
668  this%allocated_node_count = this%allocated_node_count + 1
669 
670  end associate
671  end function aabb_tree_allocate_node
672 
674  subroutine aabb_tree_deallocate_node(this, node_index)
675  class(aabb_tree_t), intent(inout) :: this
676  integer, intent(in) :: node_index
677 
678  this%nodes(node_index)%next_node_index = this%next_free_node_index
679  this%next_free_node_index = node_index
680  this%allocated_node_count = this%allocated_node_count - 1
681  end subroutine aabb_tree_deallocate_node
682 
684  subroutine aabb_tree_insert_leaf(this, leaf_node_index)
685  class(aabb_tree_t), intent(inout) :: this
686  integer, intent(in) :: leaf_node_index
687 
688  integer :: tree_node_index
689 
690  real(kind=rp) :: cost_left
691  real(kind=rp) :: cost_right
692 
693  type(aabb_node_t) :: leaf_node
694  type(aabb_node_t) :: tree_node
695  type(aabb_node_t) :: left_node
696  type(aabb_node_t) :: right_node
697 
698  type(aabb_t) :: combined_aabb
699  real(kind=rp) :: new_parent_node_cost
700  real(kind=rp) :: minimum_push_down_cost
701  type(aabb_t) :: new_left_aabb
702  type(aabb_t) :: new_right_aabb
703  integer :: leaf_sibling_index
704  type(aabb_node_t) :: leaf_sibling
705  integer :: old_parent_index
706  integer :: new_parent_index
707  type(aabb_node_t) :: new_parent
708  type(aabb_node_t) :: old_parent
709 
710  ! make sure were inserting a new leaf
711  leaf_node = this%nodes(leaf_node_index)
712 
713  ! if the tree is empty then we make the root the leaf
714  if (this%root_node_index .eq. aabb_null_node) then
715  this%root_node_index = leaf_node_index
716  leaf_node%parent_node_index = aabb_null_node
717  leaf_node%left_node_index = aabb_null_node
718  leaf_node%right_node_index = aabb_null_node
719 
720  return
721  end if
722 
723  ! search for the best place to put the new leaf in the tree
724  ! we use surface area and depth as search heuristics
725  tree_node_index = this%root_node_index
726  tree_node = this%get_node(tree_node_index)
727  do while (.not. tree_node%is_leaf())
728 
729  ! because of the test in the while loop above we know we are never a
730  ! leaf inside it
731  left_node = this%get_left_node(tree_node_index)
732  right_node = this%get_right_node(tree_node_index)
733 
734  ! ------------------------------------------------------------------- !
735 
736  combined_aabb = merge(tree_node%aabb, leaf_node%get_aabb())
737 
738  new_parent_node_cost = 2.0_rp * combined_aabb%get_surface_area()
739  minimum_push_down_cost = 2.0_rp * ( &
740  & combined_aabb%get_surface_area() &
741  & - tree_node%aabb%get_surface_area()&
742  & )
743 
744  ! use the costs to figure out whether to create a new parent here or
745  ! descend
746  if (left_node%is_leaf()) then
747  new_left_aabb = merge(leaf_node%aabb, left_node%get_aabb())
748  cost_left = new_left_aabb%get_surface_area() + minimum_push_down_cost
749  else
750  new_left_aabb = merge(leaf_node%aabb, left_node%get_aabb())
751  cost_left = ( &
752  & new_left_aabb%get_surface_area() &
753  & - left_node%aabb%get_surface_area()&
754  & ) + minimum_push_down_cost
755  end if
756 
757  if (right_node%is_leaf()) then
758  new_right_aabb = merge(leaf_node%aabb, right_node%aabb)
759  cost_right = new_right_aabb%get_surface_area() + &
760  minimum_push_down_cost
761  else
762  new_right_aabb = merge(leaf_node%aabb, right_node%aabb)
763  cost_right = ( &
764  & new_right_aabb%get_surface_area() &
765  & - right_node%aabb%get_surface_area() &
766  & ) + minimum_push_down_cost
767  end if
768 
769  ! if the cost of creating a new parent node here is less than descending
770  ! in either direction then we know we need to create a new parent node,
771  ! errrr, here and attach the leaf to that
772  if (new_parent_node_cost < cost_left .and. &
773  new_parent_node_cost < cost_right) then
774  exit
775  end if
776 
777  ! otherwise descend in the cheapest direction
778  if (cost_left .lt. cost_right) then
779  tree_node_index = tree_node%get_left_index()
780  else
781  tree_node_index = tree_node%get_right_index()
782  end if
783 
784  ! ------------------------------------------------------------------- !
785  ! Update the node and continue the loop
786  tree_node = this%get_node(tree_node_index)
787  end do
788 
789  ! the leafs sibling is going to be the node we found above and we are
790  ! going to create a new parent node and attach the leaf and this item
791  leaf_sibling_index = tree_node_index
792  leaf_sibling = this%nodes(leaf_sibling_index)
793  old_parent_index = this%get_parent_index(leaf_sibling_index)
794  new_parent_index = this%allocate_node()
795  new_parent = this%nodes(new_parent_index)
796  new_parent%parent_node_index = old_parent_index
797  new_parent%aabb = merge(leaf_node%aabb, leaf_sibling%aabb)
798 
799  if (leaf_node .lt. leaf_sibling) then
800  new_parent%left_node_index = leaf_node_index
801  new_parent%right_node_index = leaf_sibling_index
802  else
803  new_parent%left_node_index = leaf_sibling_index
804  new_parent%right_node_index = leaf_node_index
805  end if
806 
807  leaf_node%parent_node_index = new_parent_index
808  leaf_sibling%parent_node_index = new_parent_index
809 
810  if (old_parent_index .eq. aabb_null_node) then
811  ! the old parent was the root and so this is now the root
812  this%root_node_index = new_parent_index
813  else
814  ! the old parent was not the root and so we need to patch the left or
815  ! right index to point to the new node
816  old_parent = this%nodes(old_parent_index)
817  if (old_parent%left_node_index .eq. leaf_sibling_index) then
818  old_parent%left_node_index = new_parent_index
819  else
820  old_parent%right_node_index = new_parent_index
821  end if
822  this%nodes(old_parent_index) = old_parent
823  end if
824 
825  this%nodes(leaf_node_index) = leaf_node
826  this%nodes(leaf_sibling_index) = leaf_sibling
827  this%nodes(new_parent_index) = new_parent
828 
829  ! finally we need to walk back up the tree fixing heights and areas
830  tree_node_index = leaf_node%parent_node_index
831 
832  call this%fix_upwards_tree(tree_node_index)
833 
834  end subroutine aabb_tree_insert_leaf
835 
837  function aabb_tree_valid_tree(this) result(valid)
838  use stack, only: stack_i4_t
839  implicit none
840 
841  class(aabb_tree_t), intent(in) :: this
842  logical :: valid
843 
844  type(stack_i4_t) :: simple_stack
845  integer :: current_index
846  integer :: root_index, left_index, right_index
847 
848  valid = .true.
849  if (this%root_node_index .eq. aabb_null_node) then
850  valid = .false.
851  end if
852 
853  root_index = this%get_root_index()
854 
855  call simple_stack%init(this%node_capacity)
856  call simple_stack%push(root_index)
857 
858  do while (.not. simple_stack%is_empty())
859  current_index = simple_stack%pop()
860  if (current_index == aabb_null_node) cycle
861 
862  valid = valid .and. this%nodes(current_index)%is_valid()
863 
864  if (.not. this%nodes(current_index)%is_leaf()) then
865  left_index = this%get_left_index(current_index)
866  right_index = this%get_right_index(current_index)
867 
868  call simple_stack%push(left_index)
869  call simple_stack%push(right_index)
870  end if
871  end do
872  end function aabb_tree_valid_tree
873 
877  subroutine aabb_tree_fix_upwards_tree(this, tree_start_index)
878  class(aabb_tree_t), intent(inout) :: this
879  integer, intent(in) :: tree_start_index
880 
881  type(aabb_node_t) :: left_node
882  type(aabb_node_t) :: right_node
883  integer :: tree_node_index
884 
885  tree_node_index = tree_start_index
886  do while (tree_node_index .ne. aabb_null_node)
887  left_node = this%get_left_node(tree_node_index)
888  right_node = this%get_right_node(tree_node_index)
889 
890  this%nodes(tree_node_index)%aabb = merge(left_node%aabb, right_node%aabb)
891 
892  tree_node_index = this%get_parent_index(tree_node_index)
893  end do
894  end subroutine aabb_tree_fix_upwards_tree
895 
897  subroutine aabb_tree_print(this)
898  use stack, only: stack_i4_t
899  class(aabb_tree_t), intent(inout) :: this
900  type(stack_i4_t) :: simple_stack
901 
902  integer :: current_index
903  integer :: root_index, left_index, right_index
904 
905  root_index = this%get_root_index()
906  call simple_stack%init(this%node_capacity)
907  call simple_stack%push(root_index)
908 
909  do while (.not. simple_stack%is_empty())
910  current_index = simple_stack%pop()
911  if (current_index .eq. aabb_null_node) cycle
912 
913  left_index = this%get_left_index(current_index)
914  right_index = this%get_right_index(current_index)
915 
916  call simple_stack%push(this%nodes(current_index)%left_node_index)
917  call simple_stack%push(this%nodes(current_index)%right_node_index)
918 
919  write(*, *) "i = ", current_index
920  write(*, *) " Parent : ", this%get_parent_index(current_index)
921  write(*, *) " Children: ", this%get_left_index(current_index), &
922  this%get_right_index(current_index)
923 
924  write(*, *) " object_index = ", this%nodes(current_index)%object_index
925  end do
926 
927  end subroutine aabb_tree_print
928 
930  subroutine aabb_tree_resize_node_pool(this, new_capacity)
931  class(aabb_tree_t), intent(inout) :: this
932  integer, intent(in) :: new_capacity
933 
934  type(aabb_node_t), dimension(:), allocatable :: temp
935  integer :: i
936 
937  allocate(temp(new_capacity))
938  temp(:this%node_capacity) = this%nodes(:this%node_capacity)
939 
940  do i = this%allocated_node_count, new_capacity
941  temp(i)%next_node_index = i + 1
942  end do
943  temp(new_capacity)%next_node_index = aabb_null_node
944 
945  call move_alloc(temp, this%nodes)
946 
947  this%node_capacity = new_capacity
948  this%next_free_node_index = this%allocated_node_count + 1
949 
950  end subroutine aabb_tree_resize_node_pool
951 
952 end module aabb_tree
Merge two aabbs.
Definition: aabb.f90:93
Axis Aligned Bounding Box (aabb) Tree data structure.
Definition: aabb_tree.f90:70
pure type(aabb_t) function aabb_tree_get_aabb(this, node_index)
Definition: aabb_tree.f90:579
subroutine aabb_tree_fix_upwards_tree(this, tree_start_index)
Fixes the tree upwards.
Definition: aabb_tree.f90:878
pure logical function aabb_node_is_valid(this)
Returns true if the node is a valid node.
Definition: aabb_tree.f90:273
subroutine aabb_node_init(this)
Initializes the AABB node.
Definition: aabb_tree.f90:197
pure type(aabb_node_t) function aabb_tree_get_right_node(this, node_index)
Returns the right node of the node at the given index.
Definition: aabb_tree.f90:571
pure type(aabb_node_t) function aabb_tree_get_parent_node(this, node_index)
Returns the parent node of the node at the given index.
Definition: aabb_tree.f90:552
pure type(aabb_node_t) function aabb_tree_get_root_node(this)
Returns the root node of the tree.
Definition: aabb_tree.f90:543
subroutine aabb_tree_insert_object(this, object, object_index)
Inserts an object into the tree.
Definition: aabb_tree.f90:590
pure type(aabb_node_t) function aabb_tree_get_node(this, node_index)
Returns the node at the given index.
Definition: aabb_tree.f90:534
subroutine aabb_tree_print(this)
Prints the tree.
Definition: aabb_tree.f90:898
pure logical function aabb_node_greater(this, other)
Returns true if the node is greater than the other node.
Definition: aabb_tree.f90:305
integer function aabb_tree_allocate_node(this)
Allocates a new node in the tree.
Definition: aabb_tree.f90:651
pure type(aabb_t) function aabb_node_get_aabb(this)
Returns the Axis Aligned Bounding Box (aabb) of the node.
Definition: aabb_tree.f90:211
pure integer function aabb_tree_get_root_index(this)
Returns the index of the root node.
Definition: aabb_tree.f90:493
pure integer function aabb_node_get_right_index(this)
Returns the right index of the node.
Definition: aabb_tree.f90:243
subroutine aabb_tree_deallocate_node(this, node_index)
Deallocates a node in the tree.
Definition: aabb_tree.f90:675
pure integer function aabb_node_get_left_index(this)
Returns the left index of the node.
Definition: aabb_tree.f90:235
pure logical function aabb_node_is_leaf(this)
Returns true if the node is a leaf node.
Definition: aabb_tree.f90:264
real(kind=dp) function aabb_node_min_distance(this, p)
Get the minimum possible distance from the aabb to a point.
Definition: aabb_tree.f90:251
subroutine aabb_tree_init(this, initial_capacity)
Initializes the AABB tree.
Definition: aabb_tree.f90:319
subroutine aabb_tree_insert_leaf(this, leaf_node_index)
Inserts a leaf into the tree.
Definition: aabb_tree.f90:685
subroutine aabb_tree_query_overlaps(this, object, object_index, overlaps)
Queries the tree for overlapping objects.
Definition: aabb_tree.f90:605
subroutine aabb_tree_build_tree(this, objects)
Builds the tree.
Definition: aabb_tree.f90:341
integer function, dimension(:), allocatable sort(array)
Definition: aabb_tree.f90:426
pure logical function aabb_node_less(this, other)
Returns true if the node is less than the other node.
Definition: aabb_tree.f90:295
pure integer function aabb_tree_get_parent_index(this, node_index)
Returns the index of the parent node of the node at the given index.
Definition: aabb_tree.f90:502
pure integer function aabb_node_get_parent_index(this)
Returns the parent index of the node.
Definition: aabb_tree.f90:227
pure integer function aabb_tree_get_left_index(this, node_index)
Returns the index of the left node of the node at the given index.
Definition: aabb_tree.f90:512
integer, parameter, public aabb_null_node
Definition: aabb_tree.f90:78
pure integer function aabb_node_get_object_index(this)
Returns the object index of the node.
Definition: aabb_tree.f90:219
logical function aabb_tree_valid_tree(this)
Validates the tree.
Definition: aabb_tree.f90:838
pure type(aabb_node_t) function aabb_tree_get_left_node(this, node_index)
Returns the left node of the node at the given index.
Definition: aabb_tree.f90:561
integer function aabb_tree_get_size(this)
Returns the size of the tree, in number of leaves.
Definition: aabb_tree.f90:457
subroutine aabb_tree_resize_node_pool(this, new_capacity)
Resizes the node pool.
Definition: aabb_tree.f90:931
pure integer function aabb_tree_get_right_index(this, node_index)
Returns the index of the right node of the node at the given index.
Definition: aabb_tree.f90:522
Axis Aligned Bounding Box (aabb) implementation in Fortran.
Definition: aabb.f90:71
type(aabb_t) function, public get_aabb(object, padding)
Construct the aabb of a predefined object.
Definition: aabb.f90:181
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Implements a dynamic stack ADT.
Definition: stack.f90:35
Defines a triangular element.
Definition: tri.f90:34
Utilities.
Definition: utils.f90:35
Axis Aligned Bounding Box (aabb) data structure.
Definition: aabb.f90:107
Node type for the Axis Aligned Bounding Box (aabb) Tree.
Definition: aabb_tree.f90:85
Axis Aligned Bounding Box (aabb) Tree.
Definition: aabb_tree.f90:129
Integer based stack.
Definition: stack.f90:63
Triangular element.
Definition: tri.f90:60