Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
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
71 use aabb, only : aabb_t, get_aabb, merge
72 use tri, only : tri_t
73 use num_types, only : rp, dp
74 use stack, only : stack_i4_t
75 use utils, only : neko_error
76 implicit none
77 private
78
79 integer, parameter, public :: aabb_null_node = -1
80
81 ! ========================================================================== !
82 ! Type definitions
83 ! ========================================================================== !
84
86 type, public :: aabb_node_t
87 private
88 type(aabb_t), public :: aabb
89 integer :: object_index = -1
90
91 ! tree links
92 integer :: parent_node_index = aabb_null_node
93 integer :: left_node_index = aabb_null_node
94 integer :: right_node_index = aabb_null_node
95
96 ! node linked list link
97 integer :: next_node_index = aabb_null_node
98
99 contains
100 procedure, pass(this), public :: init => aabb_node_init
101
102 ! Getters
103 procedure, pass(this), public :: get_aabb => aabb_node_get_aabb
104 procedure, pass(this), public :: get_object_index => &
106 procedure, pass(this), public :: get_parent_index => &
108 procedure, pass(this), public :: get_left_index => &
110 procedure, pass(this), public :: get_right_index => &
112
113 ! Unary operations
114 procedure, pass(this), public :: min_distance => aabb_node_min_distance
115
116 ! Boolean operators
117 procedure, pass(this), public :: is_leaf => aabb_node_is_leaf
118 procedure, pass(this), public :: is_valid => aabb_node_is_valid
119
120 ! Comparison operators
121 generic :: operator(.lt.) => less
122 generic :: operator(.gt.) => greater
123
124 procedure, pass(this) :: less => aabb_node_less
125 procedure, pass(this) :: greater => aabb_node_greater
126
127 end type aabb_node_t
128
130 type, public :: aabb_tree_t
131 private
132 type(aabb_node_t), allocatable :: nodes(:)
133 integer :: root_node_index = aabb_null_node
134 integer :: allocated_node_count = 0
135 integer :: next_free_node_index = aabb_null_node
136 integer :: node_capacity = 0
137 integer :: growth_size = 1
138
139 contains
140
141 ! Initializers
142 procedure, pass(this), public :: init => aabb_tree_init
143 procedure, pass(this), public :: build => aabb_tree_build_tree
144 procedure, pass(this), public :: insert_object => &
146
147 ! Getters
148 procedure, pass(this), public :: get_size => aabb_tree_get_size
149
150 procedure, pass(this), public :: get_root_index => &
152 procedure, pass(this), public :: get_parent_index => &
154 procedure, pass(this), public :: get_left_index => &
156 procedure, pass(this), public :: get_right_index => &
158
159 procedure, pass(this), public :: get_node => aabb_tree_get_node
160 procedure, pass(this), public :: get_root_node => &
162 procedure, pass(this), public :: get_parent_node => &
164 procedure, pass(this), public :: get_left_node => &
166 procedure, pass(this), public :: get_right_node => &
168
169 procedure, pass(this), public :: get_aabb => aabb_tree_get_aabb
170
171 procedure, pass(this), public :: query_overlaps => &
173
174 procedure, pass(this), public :: print => aabb_tree_print
175
176 ! ----------------------------------------------------------------------- !
177 ! Internal methods
178
179 procedure, pass(this) :: allocate_node => aabb_tree_allocate_node
180 procedure, pass(this) :: deallocate_node => aabb_tree_deallocate_node
181 procedure, pass(this) :: resize_node_pool => aabb_tree_resize_node_pool
182 procedure, pass(this) :: insert_leaf => aabb_tree_insert_leaf
183
184 procedure, pass(this) :: fix_upwards_tree => aabb_tree_fix_upwards_tree
185
186 procedure, pass(this) :: valid_tree => aabb_tree_valid_tree
187
188 end type aabb_tree_t
189
190contains
191
192 ! ========================================================================== !
193 ! Definitions of node methods
194 ! ========================================================================== !
195
197 subroutine aabb_node_init(this)
198 class(aabb_node_t), intent(inout) :: this
199
200 this%object_index = -1
201 this%parent_node_index = aabb_null_node
202 this%left_node_index = aabb_null_node
203 this%right_node_index = aabb_null_node
204 this%next_node_index = aabb_null_node
205 end subroutine aabb_node_init
206
207 ! -------------------------------------------------------------------------- !
208 ! Getters
209
211 pure function aabb_node_get_aabb(this) result(res)
212 class(aabb_node_t), intent(in) :: this
213 type(aabb_t) :: res
214
215 res = this%aabb
216 end function aabb_node_get_aabb
217
219 pure function aabb_node_get_object_index(this) result(object_index)
220 class(aabb_node_t), intent(in) :: this
221 integer :: object_index
222
223 object_index = this%object_index
224 end function aabb_node_get_object_index
225
227 pure function aabb_node_get_parent_index(this) result(parent_index)
228 class(aabb_node_t), intent(in) :: this
229 integer :: parent_index
230
231 parent_index = this%parent_node_index
232 end function aabb_node_get_parent_index
233
235 pure function aabb_node_get_left_index(this) result(left_index)
236 class(aabb_node_t), intent(in) :: this
237 integer :: left_index
238
239 left_index = this%left_node_index
240 end function aabb_node_get_left_index
241
243 pure function aabb_node_get_right_index(this) result(right_index)
244 class(aabb_node_t), intent(in) :: this
245 integer :: right_index
246
247 right_index = this%right_node_index
248 end function aabb_node_get_right_index
249
251 function aabb_node_min_distance(this, p) result(distance)
252 class(aabb_node_t), intent(in) :: this
253 real(kind=dp), dimension(3), intent(in) :: p
254 real(kind=dp) :: distance
255
256 distance = 0.5_rp * this%aabb%get_diameter() &
257 - norm2(this%aabb%get_center() - p)
258 end function aabb_node_min_distance
259
260 ! -------------------------------------------------------------------------- !
261 ! Boolean operators
262
264 pure function aabb_node_is_leaf(this) result(res)
265 class(aabb_node_t), intent(in) :: this
266 logical :: res
267
268 res = this%left_node_index == aabb_null_node .and. &
269 this%right_node_index == aabb_null_node
270 end function aabb_node_is_leaf
271
273 pure function aabb_node_is_valid(this) result(valid)
274 class(aabb_node_t), intent(in) :: this
275 logical :: valid
276
277 if (this%is_leaf()) then
278 valid = &
279 & this%left_node_index .eq. aabb_null_node .and. &
280 & this%right_node_index .eq. aabb_null_node .and. &
281 & this%object_index .gt. 0
282 else
283 valid = &
284 & this%left_node_index .ne. aabb_null_node .and. &
285 & this%right_node_index .ne. aabb_null_node .and. &
286 & this%object_index .eq. -1
287 end if
288
289 end function aabb_node_is_valid
290
291 ! -------------------------------------------------------------------------- !
292 ! Comparison operators
293
295 pure function aabb_node_less(this, other) result(res)
296 class(aabb_node_t), intent(in) :: this
297 class(aabb_node_t), intent(in) :: other
298 logical :: res
299
300 res = this%aabb .lt. other%aabb
301
302 end function aabb_node_less
303
305 pure function aabb_node_greater(this, other) result(res)
306 class(aabb_node_t), intent(in) :: this
307 class(aabb_node_t), intent(in) :: other
308 logical :: res
309
310 res = this%aabb .gt. other%aabb
311
312 end function aabb_node_greater
313
314 ! ========================================================================== !
315 ! Definitions of tree methods
316 ! ========================================================================== !
317
319 subroutine aabb_tree_init(this, initial_capacity)
320 class(aabb_tree_t), intent(inout) :: this
321 integer, intent(in) :: initial_capacity
322
323 integer :: i
324
325 this%root_node_index = aabb_null_node
326 this%allocated_node_count = 0
327 this%next_free_node_index = 1
328 this%node_capacity = initial_capacity
329 this%growth_size = initial_capacity
330
331 if (allocated(this%nodes)) deallocate(this%nodes)
332 allocate(this%nodes(initial_capacity))
333
334 do i = 1, initial_capacity
335 this%nodes(i)%next_node_index = i + 1
336 end do
337 this%nodes(initial_capacity)%next_node_index = aabb_null_node
338 end subroutine aabb_tree_init
339
341 subroutine aabb_tree_build_tree(this, objects)
342 class(aabb_tree_t), intent(inout) :: this
343 class(*), dimension(:), intent(in) :: objects
344
345 integer :: i_obj, i_node, i
346 logical :: done
347
348 integer :: start_layer, end_layer
349
350 type(aabb_t), dimension(:), allocatable :: box_list
351 integer, dimension(:), allocatable :: sorted_indices
352
353 call this%init(size(objects) * 2)
354
355 ! ------------------------------------------------------------------------ !
356 ! Start by sorting the list of objects, then build a balanced binary tree
357 ! from the sorted list
358
359 allocate(box_list(size(objects)))
360
361 do i_obj = 1, size(objects)
362 box_list(i_obj) = get_aabb(objects(i_obj))
363 end do
364 call sort(box_list, sorted_indices)
365
366 do i = 1, size(sorted_indices)
367 i_obj = sorted_indices(i)
368 i_node = this%allocate_node()
369 this%nodes(i_node)%aabb = get_aabb(objects(i_obj))
370 this%nodes(i_node)%object_index = i_obj
371 end do
372
373
374 start_layer = 1
375 end_layer = size(objects)
376 done = .false.
377 do while (.not. done)
378
379 ! build the next layer
380 do i = start_layer, end_layer - 1, 2
381 i_node = this%allocate_node()
382
383 this%nodes(i_node)%aabb = merge(this%nodes(i)%aabb, &
384 this%nodes(i + 1)%aabb)
385
386 this%nodes(i_node)%left_node_index = i
387 this%nodes(i_node)%right_node_index = i + 1
388
389 this%nodes(i)%parent_node_index = i_node
390 this%nodes(i + 1)%parent_node_index = i_node
391 end do
392
393 ! if the number of nodes is odd, we need to create a new node to hold the
394 ! last node
395 if (mod(end_layer - start_layer, 2) .eq. 0) then
396 i_node = this%allocate_node()
397 this%nodes(i_node)%aabb = this%nodes(end_layer)%aabb
398 this%nodes(i_node)%left_node_index = end_layer
399 this%nodes(i_node)%right_node_index = aabb_null_node
400
401 this%nodes(end_layer)%parent_node_index = i_node
402 end if
403
404 ! move to the next layer
405 start_layer = end_layer + 1
406 end_layer = this%allocated_node_count
407
408 ! If there is only one node left, we are done
409 done = start_layer .eq. end_layer
410 end do
411
412 ! The last node allocated is the root node
413 this%root_node_index = this%allocated_node_count
414
415 if (this%get_size() .ne. size(objects)) then
416 print *, "this%get_size() = ", this%get_size()
417 print *, "size(objects) = ", size(objects)
418 call neko_error("Invalid tree size")
419 end if
420
421 end subroutine aabb_tree_build_tree
422
424 subroutine sort(array, indices)
425 type(aabb_t), dimension(:), intent(in) :: array
426 integer, intent(inout), dimension(:), allocatable :: indices
427 logical, dimension(:), allocatable :: visited
428
429 integer :: i, imin
430 integer :: minidx
431
432 allocate(indices(size(array)))
433 allocate(visited(size(array)))
434
435 visited = .false.
436 indices = 0
437 do i = 1, size(array)
438 minidx = -1
439 do imin = 1, size(array)
440 if (.not. visited(imin) .and. minidx .eq. -1) minidx = imin
441 if (minidx .gt. -1) then
442 if (visited(imin) .and. array(imin) .lt. array(minidx)) minidx = imin
443 end if
444 end do
445
446 indices(i) = minidx
447 visited(minidx) = .true.
448 end do
449
450 end subroutine sort
451
452 ! -------------------------------------------------------------------------- !
453 ! Getters
454
456 function aabb_tree_get_size(this) result(size)
457 class(aabb_tree_t), intent(in) :: this
458 integer :: size
459
460 type(stack_i4_t) :: simple_stack
461 integer :: idx, tmp
462
463 call simple_stack%init(this%allocated_node_count)
464 size = 0
465 tmp = this%get_root_index()
466 if (tmp .ne. aabb_null_node) then
467 call simple_stack%push(tmp)
468 end if
469
470 do while (.not. simple_stack%is_empty())
471 idx = simple_stack%pop()
472 if (idx .eq. aabb_null_node) cycle
473
474 if (this%nodes(idx)%is_leaf()) then
475 size = size + 1
476 else
477 tmp = this%get_left_index(idx)
478 call simple_stack%push(tmp)
479 tmp = this%get_right_index(idx)
480 call simple_stack%push(tmp)
481 end if
482 end do
483
484 end function aabb_tree_get_size
485
486 ! -------------------------------------------------------------------------- !
487 ! Get index of nodes
488
490 pure function aabb_tree_get_root_index(this) result(root_index)
491 class(aabb_tree_t), intent(in) :: this
492 integer :: root_index
493
494 root_index = this%root_node_index
495 end function aabb_tree_get_root_index
496
498 pure function aabb_tree_get_parent_index(this, node_index) &
499 result(parent_index)
500 class(aabb_tree_t), intent(in) :: this
501 integer, intent(in) :: node_index
502 integer :: parent_index
503
504 parent_index = this%nodes(node_index)%parent_node_index
505 end function aabb_tree_get_parent_index
506
508 pure function aabb_tree_get_left_index(this, node_index) &
509 result(left_index)
510 class(aabb_tree_t), intent(in) :: this
511 integer, intent(in) :: node_index
512 integer :: left_index
513
514 left_index = this%nodes(node_index)%left_node_index
515 end function aabb_tree_get_left_index
516
518 pure function aabb_tree_get_right_index(this, node_index) &
519 result(right_index)
520 class(aabb_tree_t), intent(in) :: this
521 integer, intent(in) :: node_index
522 integer :: right_index
523
524 right_index = this%nodes(node_index)%right_node_index
525 end function aabb_tree_get_right_index
526
527 ! -------------------------------------------------------------------------- !
528 ! Get nodes
529
531 pure function aabb_tree_get_node(this, node_index) result(node)
532 class(aabb_tree_t), intent(in) :: this
533 integer, intent(in) :: node_index
534 type(aabb_node_t) :: node
535
536 node = this%nodes(node_index)
537 end function aabb_tree_get_node
538
540 pure function aabb_tree_get_root_node(this) result(root_node)
541 class(aabb_tree_t), intent(in) :: this
542 type(aabb_node_t) :: root_node
543
544 root_node = this%nodes(this%root_node_index)
545 end function aabb_tree_get_root_node
546
548 pure function aabb_tree_get_parent_node(this, node_index) &
549 result(parent_node)
550 class(aabb_tree_t), intent(in) :: this
551 integer, intent(in) :: node_index
552 type(aabb_node_t) :: parent_node
553
554 parent_node = this%nodes(this%nodes(node_index)%parent_node_index)
555 end function aabb_tree_get_parent_node
556
558 pure function aabb_tree_get_left_node(this, node_index) result(left_node)
559 class(aabb_tree_t), intent(in) :: this
560 integer, intent(in) :: node_index
561 type(aabb_node_t) :: left_node
562
563 left_node = this%nodes(this%nodes(node_index)%left_node_index)
564 end function aabb_tree_get_left_node
565
567 pure function aabb_tree_get_right_node(this, node_index) &
568 result(right_node)
569 class(aabb_tree_t), intent(in) :: this
570 integer, intent(in) :: node_index
571 type(aabb_node_t) :: right_node
572
573 right_node = this%nodes(this%nodes(node_index)%right_node_index)
574 end function aabb_tree_get_right_node
575
576 pure function aabb_tree_get_aabb(this, node_index) result(out_box)
577 class(aabb_tree_t), intent(in) :: this
578 integer, intent(in) :: node_index
579 type(aabb_t) :: out_box
580
581 out_box = this%nodes(node_index)%aabb
582 end function aabb_tree_get_aabb
583
584 ! -------------------------------------------------------------------------- !
585
587 subroutine aabb_tree_insert_object(this, object, object_index)
588 class(aabb_tree_t), intent(inout) :: this
589 class(*), intent(in) :: object
590 integer, intent(in) :: object_index
591
592 integer :: node_index
593
594 node_index = this%allocate_node()
595 this%nodes(node_index)%aabb = get_aabb(object)
596 this%nodes(node_index)%object_index = object_index
597
598 call this%insert_leaf(node_index)
599 end subroutine aabb_tree_insert_object
600
602 subroutine aabb_tree_query_overlaps(this, object, object_index, overlaps)
603 class(aabb_tree_t), intent(in) :: this
604 class(*), intent(in) :: object
605 integer, intent(in) :: object_index
606 integer, intent(inout), allocatable :: overlaps(:)
607
608 type(stack_i4_t) :: simple_stack
609 type(aabb_t) :: object_box
610
611 integer :: root_index, left_index, right_index
612 integer :: node_index
613
614 if (allocated(overlaps)) deallocate(overlaps)
615 allocate(overlaps(0))
616
617 object_box = get_aabb(object)
618 root_index = this%get_root_index()
619
620 call simple_stack%init()
621 call simple_stack%push(root_index)
622
623 do while (.not. simple_stack%is_empty())
624 node_index = simple_stack%pop()
625
626 if (node_index == aabb_null_node) cycle
627
628 if (this%nodes(node_index)%aabb%overlaps(object_box)) then
629 if (this%nodes(node_index)%is_leaf()) then
630 if (this%nodes(node_index)%object_index .ne. object_index) then
631 overlaps = [this%nodes(node_index)%object_index, overlaps]
632 end if
633 else
634 left_index = this%get_left_index(node_index)
635 call simple_stack%push(left_index)
636 right_index = this%get_right_index(node_index)
637 call simple_stack%push(right_index)
638 end if
639 end if
640 end do
641 call simple_stack%free()
642 end subroutine aabb_tree_query_overlaps
643
644 ! -------------------------------------------------------------------------- !
645 ! Internal methods
646
648 function aabb_tree_allocate_node(this) result(node_index)
649 class(aabb_tree_t), intent(inout) :: this
650 integer :: node_index
651
652 if (this%next_free_node_index == aabb_null_node) then
653 call this%resize_node_pool(this%node_capacity + this%growth_size)
654 end if
655
656 node_index = this%next_free_node_index
657
658 associate(new_node => this%nodes(node_index))
659 this%next_free_node_index = new_node%next_node_index
660
661 new_node%parent_node_index = aabb_null_node
662 new_node%left_node_index = aabb_null_node
663 new_node%right_node_index = aabb_null_node
664
665 this%next_free_node_index = new_node%next_node_index
666 this%allocated_node_count = this%allocated_node_count + 1
667
668 end associate
669 end function aabb_tree_allocate_node
670
672 subroutine aabb_tree_deallocate_node(this, node_index)
673 class(aabb_tree_t), intent(inout) :: this
674 integer, intent(in) :: node_index
675
676 this%nodes(node_index)%next_node_index = this%next_free_node_index
677 this%next_free_node_index = node_index
678 this%allocated_node_count = this%allocated_node_count - 1
679 end subroutine aabb_tree_deallocate_node
680
682 subroutine aabb_tree_insert_leaf(this, leaf_node_index)
683 class(aabb_tree_t), intent(inout) :: this
684 integer, intent(in) :: leaf_node_index
685
686 integer :: tree_node_index
687
688 real(kind=rp) :: cost_left
689 real(kind=rp) :: cost_right
690
691 type(aabb_node_t) :: leaf_node
692 type(aabb_node_t) :: tree_node
693 type(aabb_node_t) :: left_node
694 type(aabb_node_t) :: right_node
695
696 type(aabb_t) :: combined_aabb
697 real(kind=rp) :: new_parent_node_cost
698 real(kind=rp) :: minimum_push_down_cost
699 type(aabb_t) :: new_left_aabb
700 type(aabb_t) :: new_right_aabb
701 integer :: leaf_sibling_index
702 type(aabb_node_t) :: leaf_sibling
703 integer :: old_parent_index
704 integer :: new_parent_index
705 type(aabb_node_t) :: new_parent
706 type(aabb_node_t) :: old_parent
707
708 ! make sure were inserting a new leaf
709 leaf_node = this%nodes(leaf_node_index)
710
711 ! if the tree is empty then we make the root the leaf
712 if (this%root_node_index .eq. aabb_null_node) then
713 this%root_node_index = leaf_node_index
714 leaf_node%parent_node_index = aabb_null_node
715 leaf_node%left_node_index = aabb_null_node
716 leaf_node%right_node_index = aabb_null_node
717
718 return
719 end if
720
721 ! search for the best place to put the new leaf in the tree
722 ! we use surface area and depth as search heuristics
723 tree_node_index = this%root_node_index
724 tree_node = this%get_node(tree_node_index)
725 do while (.not. tree_node%is_leaf())
726
727 ! because of the test in the while loop above we know we are never a
728 ! leaf inside it
729 left_node = this%get_left_node(tree_node_index)
730 right_node = this%get_right_node(tree_node_index)
731
732 ! ------------------------------------------------------------------- !
733
734 combined_aabb = merge(tree_node%aabb, leaf_node%get_aabb())
735
736 new_parent_node_cost = 2.0_rp * combined_aabb%get_surface_area()
737 minimum_push_down_cost = 2.0_rp * ( &
738 & combined_aabb%get_surface_area() &
739 & - tree_node%aabb%get_surface_area()&
740 & )
741
742 ! use the costs to figure out whether to create a new parent here or
743 ! descend
744 if (left_node%is_leaf()) then
745 new_left_aabb = merge(leaf_node%aabb, left_node%get_aabb())
746 cost_left = new_left_aabb%get_surface_area() + minimum_push_down_cost
747 else
748 new_left_aabb = merge(leaf_node%aabb, left_node%get_aabb())
749 cost_left = ( &
750 & new_left_aabb%get_surface_area() &
751 & - left_node%aabb%get_surface_area()&
752 & ) + minimum_push_down_cost
753 end if
754
755 if (right_node%is_leaf()) then
756 new_right_aabb = merge(leaf_node%aabb, right_node%aabb)
757 cost_right = new_right_aabb%get_surface_area() + &
758 minimum_push_down_cost
759 else
760 new_right_aabb = merge(leaf_node%aabb, right_node%aabb)
761 cost_right = ( &
762 & new_right_aabb%get_surface_area() &
763 & - right_node%aabb%get_surface_area() &
764 & ) + minimum_push_down_cost
765 end if
766
767 ! if the cost of creating a new parent node here is less than descending
768 ! in either direction then we know we need to create a new parent node,
769 ! errrr, here and attach the leaf to that
770 if (new_parent_node_cost < cost_left .and. &
771 new_parent_node_cost < cost_right) then
772 exit
773 end if
774
775 ! otherwise descend in the cheapest direction
776 if (cost_left .lt. cost_right) then
777 tree_node_index = tree_node%get_left_index()
778 else
779 tree_node_index = tree_node%get_right_index()
780 end if
781
782 ! ------------------------------------------------------------------- !
783 ! Update the node and continue the loop
784 tree_node = this%get_node(tree_node_index)
785 end do
786
787 ! the leafs sibling is going to be the node we found above and we are
788 ! going to create a new parent node and attach the leaf and this item
789 leaf_sibling_index = tree_node_index
790 leaf_sibling = this%nodes(leaf_sibling_index)
791 old_parent_index = this%get_parent_index(leaf_sibling_index)
792 new_parent_index = this%allocate_node()
793 new_parent = this%nodes(new_parent_index)
794 new_parent%parent_node_index = old_parent_index
795 new_parent%aabb = merge(leaf_node%aabb, leaf_sibling%aabb)
796
797 if (leaf_node .lt. leaf_sibling) then
798 new_parent%left_node_index = leaf_node_index
799 new_parent%right_node_index = leaf_sibling_index
800 else
801 new_parent%left_node_index = leaf_sibling_index
802 new_parent%right_node_index = leaf_node_index
803 end if
804
805 leaf_node%parent_node_index = new_parent_index
806 leaf_sibling%parent_node_index = new_parent_index
807
808 if (old_parent_index .eq. aabb_null_node) then
809 ! the old parent was the root and so this is now the root
810 this%root_node_index = new_parent_index
811 else
812 ! the old parent was not the root and so we need to patch the left or
813 ! right index to point to the new node
814 old_parent = this%nodes(old_parent_index)
815 if (old_parent%left_node_index .eq. leaf_sibling_index) then
816 old_parent%left_node_index = new_parent_index
817 else
818 old_parent%right_node_index = new_parent_index
819 end if
820 this%nodes(old_parent_index) = old_parent
821 end if
822
823 this%nodes(leaf_node_index) = leaf_node
824 this%nodes(leaf_sibling_index) = leaf_sibling
825 this%nodes(new_parent_index) = new_parent
826
827 ! finally we need to walk back up the tree fixing heights and areas
828 tree_node_index = leaf_node%parent_node_index
829
830 call this%fix_upwards_tree(tree_node_index)
831
832 end subroutine aabb_tree_insert_leaf
833
835 function aabb_tree_valid_tree(this) result(valid)
836 class(aabb_tree_t), intent(in) :: this
837 logical :: valid
838
839 type(stack_i4_t) :: simple_stack
840 integer :: current_index
841 integer :: root_index, left_index, right_index
842
843 valid = .true.
844 if (this%root_node_index .eq. aabb_null_node) then
845 valid = .false.
846 end if
847
848 root_index = this%get_root_index()
849
850 call simple_stack%init(this%node_capacity)
851 call simple_stack%push(root_index)
852
853 do while (.not. simple_stack%is_empty())
854 current_index = simple_stack%pop()
855 if (current_index == aabb_null_node) cycle
856
857 valid = valid .and. this%nodes(current_index)%is_valid()
858
859 if (.not. this%nodes(current_index)%is_leaf()) then
860 left_index = this%get_left_index(current_index)
861 right_index = this%get_right_index(current_index)
862
863 call simple_stack%push(left_index)
864 call simple_stack%push(right_index)
865 end if
866 end do
867 end function aabb_tree_valid_tree
868
872 subroutine aabb_tree_fix_upwards_tree(this, tree_start_index)
873 class(aabb_tree_t), intent(inout) :: this
874 integer, intent(in) :: tree_start_index
875
876 type(aabb_node_t) :: left_node
877 type(aabb_node_t) :: right_node
878 integer :: tree_node_index
879
880 tree_node_index = tree_start_index
881 do while (tree_node_index .ne. aabb_null_node)
882 left_node = this%get_left_node(tree_node_index)
883 right_node = this%get_right_node(tree_node_index)
884
885 this%nodes(tree_node_index)%aabb = merge(left_node%aabb, right_node%aabb)
886
887 tree_node_index = this%get_parent_index(tree_node_index)
888 end do
889 end subroutine aabb_tree_fix_upwards_tree
890
892 subroutine aabb_tree_print(this)
893 class(aabb_tree_t), intent(inout) :: this
894 type(stack_i4_t) :: simple_stack
895
896 integer :: current_index
897 integer :: root_index, left_index, right_index
898
899 root_index = this%get_root_index()
900 call simple_stack%init(this%node_capacity)
901 call simple_stack%push(root_index)
902
903 do while (.not. simple_stack%is_empty())
904 current_index = simple_stack%pop()
905 if (current_index .eq. aabb_null_node) cycle
906
907 left_index = this%get_left_index(current_index)
908 right_index = this%get_right_index(current_index)
909
910 call simple_stack%push(this%nodes(current_index)%left_node_index)
911 call simple_stack%push(this%nodes(current_index)%right_node_index)
912
913 write(*, *) "i = ", current_index
914 write(*, *) " Parent : ", this%get_parent_index(current_index)
915 write(*, *) " Children: ", this%get_left_index(current_index), &
916 this%get_right_index(current_index)
917
918 write(*, *) " object_index = ", this%nodes(current_index)%object_index
919 end do
920
921 end subroutine aabb_tree_print
922
924 subroutine aabb_tree_resize_node_pool(this, new_capacity)
925 class(aabb_tree_t), intent(inout) :: this
926 integer, intent(in) :: new_capacity
927
928 type(aabb_node_t), dimension(:), allocatable :: temp
929 integer :: i
930
931 allocate(temp(new_capacity))
932 temp(:this%node_capacity) = this%nodes(:this%node_capacity)
933
934 do i = this%allocated_node_count, new_capacity
935 temp(i)%next_node_index = i + 1
936 end do
937 temp(new_capacity)%next_node_index = aabb_null_node
938
939 call move_alloc(temp, this%nodes)
940
941 this%node_capacity = new_capacity
942 this%next_free_node_index = this%allocated_node_count + 1
943
944 end subroutine aabb_tree_resize_node_pool
945
946end 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)
subroutine aabb_tree_fix_upwards_tree(this, tree_start_index)
Fixes the tree upwards.
pure logical function aabb_node_is_valid(this)
Returns true if the node is a valid node.
subroutine aabb_node_init(this)
Initializes the AABB node.
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.
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.
pure type(aabb_node_t) function aabb_tree_get_root_node(this)
Returns the root node of the tree.
subroutine aabb_tree_insert_object(this, object, object_index)
Inserts an object into the tree.
pure type(aabb_node_t) function aabb_tree_get_node(this, node_index)
Returns the node at the given index.
subroutine aabb_tree_print(this)
Prints the tree.
pure logical function aabb_node_greater(this, other)
Returns true if the node is greater than the other node.
integer function aabb_tree_allocate_node(this)
Allocates a new node in the tree.
pure type(aabb_t) function aabb_node_get_aabb(this)
Returns the Axis Aligned Bounding Box (aabb) of the node.
pure integer function aabb_tree_get_root_index(this)
Returns the index of the root node.
pure integer function aabb_node_get_right_index(this)
Returns the right index of the node.
subroutine aabb_tree_deallocate_node(this, node_index)
Deallocates a node in the tree.
pure integer function aabb_node_get_left_index(this)
Returns the left index of the node.
pure logical function aabb_node_is_leaf(this)
Returns true if the node is a leaf node.
real(kind=dp) function aabb_node_min_distance(this, p)
Get the minimum possible distance from the aabb to a point.
subroutine aabb_tree_init(this, initial_capacity)
Initializes the AABB tree.
subroutine aabb_tree_insert_leaf(this, leaf_node_index)
Inserts a leaf into the tree.
subroutine aabb_tree_query_overlaps(this, object, object_index, overlaps)
Queries the tree for overlapping objects.
subroutine aabb_tree_build_tree(this, objects)
Builds the tree.
subroutine sort(array, indices)
Return a list of sorted indices of the aabb nodes.
pure logical function aabb_node_less(this, other)
Returns true if the node is less than the other node.
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.
pure integer function aabb_node_get_parent_index(this)
Returns the parent index of the node.
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.
integer, parameter, public aabb_null_node
Definition aabb_tree.f90:79
pure integer function aabb_node_get_object_index(this)
Returns the object index of the node.
logical function aabb_tree_valid_tree(this)
Validates the tree.
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.
integer function aabb_tree_get_size(this)
Returns the size of the tree, in number of leaves.
subroutine aabb_tree_resize_node_pool(this, new_capacity)
Resizes the node pool.
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.
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:86
Axis Aligned Bounding Box (aabb) Tree.
Integer based stack.
Definition stack.f90:63
Triangular element.
Definition tri.f90:61