Neko  0.9.0
A portable framework for high-order spectral element flow simulations
aabb.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/include/AABB.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 module aabb
72  use num_types, only: dp
73  use element, only: element_t
74  use point, only: point_t
75  use tri, only: tri_t
76  use quad, only: quad_t
77  use tet, only: tet_t
78  use hex, only: hex_t
79  use mesh, only: mesh_t
80  use tri_mesh, only: tri_mesh_t
81  use tet_mesh, only: tet_mesh_t
82  use utils, only: neko_error
83 
84  implicit none
85  private
86  public :: aabb_t, get_aabb, merge, intersection
87 
88  ! ========================================================================== !
89  ! Public interface for free functions
90  ! ========================================================================== !
91 
93  interface merge
94  module procedure merge_aabb
95  end interface merge
96 
98  interface intersection
99  module procedure intersection_aabb
100  end interface intersection
101 
107  type :: aabb_t
108  private
109 
110  logical :: initialized = .false.
111  real(kind=dp) :: box_min(3) = huge(0.0_dp)
112  real(kind=dp) :: box_max(3) = -huge(0.0_dp)
113  real(kind=dp) :: center(3) = 0.0_dp
114  real(kind=dp) :: diameter = huge(0.0_dp)
115  real(kind=dp) :: surface_area = 0.0_dp
116 
117  contains
118 
119  ! Initializers
120  procedure, pass(this), public :: init => aabb_init
121 
122  ! Getters
123  procedure, pass(this), public :: get_min => aabb_get_min
124  procedure, pass(this), public :: get_max => aabb_get_max
125  procedure, pass(this), public :: get_width => aabb_get_width
126  procedure, pass(this), public :: get_height => aabb_get_height
127  procedure, pass(this), public :: get_depth => aabb_get_depth
128  procedure, pass(this), public :: get_diameter => aabb_get_diameter
129  procedure, pass(this), public :: get_surface_area => aabb_get_surface_area
130  procedure, pass(this), public :: get_center => aabb_get_center
131  procedure, pass(this), public :: get_diagonal => aabb_get_diagonal
132 
133  procedure, pass(this), public :: add_padding
134 
135  ! Comparison operators
136  generic :: operator(.lt.) => less
137  generic :: operator(.gt.) => greater
138 
140  procedure, pass(this), public :: overlaps => aabb_overlaps
142  procedure, pass(this), public :: contains => aabb_contains_other
144  procedure, pass(this), public :: contains_point => aabb_contains_point
145 
146  ! Private comparison operators
147  procedure, pass(this) :: less => aabb_less
148  procedure, pass(this) :: greater => aabb_greater
149 
150  ! Private operations
151  procedure, pass(this), private :: calculate_surface_area
152 
153  end type aabb_t
154 
155 contains
156 
157  ! ========================================================================== !
158  ! Constructors
159  ! ========================================================================== !
160 
180  function get_aabb(object, padding) result(box)
181  use utils, only: neko_error
182  implicit none
183 
184  class(*), intent(in) :: object
185  real(kind=dp), intent(in), optional :: padding
186  type(aabb_t) :: box
187 
188  select type(object)
189 
190  type is (tri_t)
191  box = get_aabb_element(object)
192  type is (hex_t)
193  box = get_aabb_element(object)
194  type is (tet_t)
195  box = get_aabb_element(object)
196  type is (quad_t)
197  box = get_aabb_element(object)
198 
199  type is (mesh_t)
200  box = get_aabb_mesh(object)
201  type is (tri_mesh_t)
202  box = get_aabb_tri_mesh(object)
203  type is (tet_mesh_t)
204  box = get_aabb_tet_mesh(object)
205 
206  class default
207  call neko_error("get_aabb: Unsupported object type")
208  end select
209 
210  if (present(padding)) then
211  call box%add_padding(padding)
212  end if
213 
214  end function get_aabb
215 
221  subroutine add_padding(this, padding)
222  class(aabb_t), intent(inout) :: this
223  real(kind=dp), intent(in) :: padding
224  real(kind=dp) :: box_min(3), box_max(3)
225 
226  box_min = this%box_min - padding * (this%diameter)
227  box_max = this%box_max + padding * (this%diameter)
228 
229  call this%init(box_min, box_max)
230  end subroutine add_padding
231 
241  function get_aabb_element(object) result(box)
242  class(element_t), intent(in) :: object
243  type(aabb_t) :: box
244 
245  integer :: i
246  type(point_t), pointer :: pi
247  real(kind=dp) :: box_min(3), box_max(3)
248 
249  box_min = huge(0.0_dp); box_max = -huge(0.0_dp)
250 
251  do i = 1, object%n_points()
252  pi => object%p(i)
253  box_min = min(box_min, pi%x)
254  box_max = max(box_max, pi%x)
255  end do
256 
257  call box%init(box_min, box_max)
258  end function get_aabb_element
259 
268  function get_aabb_mesh(object) result(box)
269  type(mesh_t), intent(in) :: object
270  type(aabb_t) :: box
271 
272  integer :: i
273  type(aabb_t) :: temp_box
274 
275  do i = 1, object%nelv
276  temp_box = get_aabb(object%elements(i))
277  box = merge(box, temp_box)
278  end do
279 
280  end function get_aabb_mesh
281 
290  function get_aabb_tri_mesh(object) result(box)
291  type(tri_mesh_t), intent(in) :: object
292  type(aabb_t) :: box
293 
294  integer :: i
295  type(aabb_t) :: temp_box
296 
297  do i = 1, object%nelv
298  temp_box = get_aabb(object%el(i))
299  box = merge(box, temp_box)
300  end do
301 
302  end function get_aabb_tri_mesh
303 
312  function get_aabb_tet_mesh(object) result(box)
313  type(tet_mesh_t), intent(in) :: object
314  type(aabb_t) :: box
315 
316  integer :: i
317  type(aabb_t) :: temp_box
318 
319  do i = 1, object%nelv
320  temp_box = get_aabb(object%el(i))
321  box = merge(box, temp_box)
322  end do
323 
324  end function get_aabb_tet_mesh
325 
326 
327  ! ========================================================================== !
328  ! Initializers
329  ! ========================================================================== !
330 
334  subroutine aabb_init(this, lower_left_front, upper_right_back)
335  class(aabb_t), intent(inout) :: this
336  real(kind=dp), dimension(3), intent(in) :: lower_left_front
337  real(kind=dp), dimension(3), intent(in) :: upper_right_back
338 
339  this%initialized = .true.
340  this%box_min = lower_left_front
341  this%box_max = upper_right_back
342  this%center = (this%box_min + this%box_max) / 2.0_dp
343  this%diameter = norm2(this%box_max - this%box_min)
344 
345  this%surface_area = this%calculate_surface_area()
346 
347  end subroutine aabb_init
348 
349  ! ========================================================================== !
350  ! Getters
351  ! ========================================================================== !
352 
354  pure function aabb_get_min(this) result(min)
355  class(aabb_t), intent(in) :: this
356  real(kind=dp), dimension(3) :: min
357 
358  min = this%box_min
359  end function aabb_get_min
360 
362  pure function aabb_get_max(this) result(max)
363  class(aabb_t), intent(in) :: this
364  real(kind=dp), dimension(3) :: max
365 
366  max = this%box_max
367  end function aabb_get_max
368 
370  pure function aabb_get_width(this) result(width)
371  class(aabb_t), intent(in) :: this
372  real(kind=dp) :: width
373 
374  width = this%box_max(1) - this%box_min(1)
375  end function aabb_get_width
376 
378  pure function aabb_get_depth(this) result(depth)
379  class(aabb_t), intent(in) :: this
380  real(kind=dp) :: depth
381 
382  depth = this%box_max(2) - this%box_min(2)
383  end function aabb_get_depth
384 
386  pure function aabb_get_height(this) result(height)
387  class(aabb_t), intent(in) :: this
388  real(kind=dp) :: height
389 
390  height = this%box_max(3) - this%box_min(3)
391  end function aabb_get_height
392 
394  pure function aabb_get_diameter(this) result(diameter)
395  class(aabb_t), intent(in) :: this
396  real(kind=dp) :: diameter
397 
398  diameter = this%diameter
399  end function aabb_get_diameter
400 
402  pure function aabb_get_surface_area(this) result(surface_area)
403  class(aabb_t), intent(in) :: this
404  real(kind=dp) :: surface_area
405 
406  surface_area = this%surface_area
407  end function aabb_get_surface_area
408 
410  pure function aabb_get_center(this) result(center)
411  class(aabb_t), intent(in) :: this
412  real(kind=dp), dimension(3) :: center
413 
414  center = this%center
415  end function aabb_get_center
416 
418  pure function aabb_get_diagonal(this) result(diagonal)
419  class(aabb_t), intent(in) :: this
420  real(kind=dp), dimension(3) :: diagonal
421 
422  diagonal = this%box_max - this%box_min
423  end function aabb_get_diagonal
424 
425  ! ========================================================================== !
426  ! Operations
427  ! ========================================================================== !
428 
430  function aabb_overlaps(this, other) result(is_overlapping)
431  class(aabb_t), intent(in) :: this
432  class(aabb_t), intent(in) :: other
433  logical :: is_overlapping
434 
435  if (.not. this%initialized .or. .not. other%initialized) then
436  ! call neko_error("aabb_overlaps: One or both aabbs are not initialized")
437  is_overlapping = .false.
438  else
439 
440  is_overlapping = all(this%box_min .le. other%box_max) .and. &
441  all(this%box_max .ge. other%box_min)
442  end if
443 
444  end function aabb_overlaps
445 
447  function aabb_contains_other(this, other) result(is_contained)
448  class(aabb_t), intent(in) :: this
449  class(aabb_t), intent(in) :: other
450  logical :: is_contained
451 
452  ! if (.not. this%initialized .or. .not. other%initialized) then
453  ! call neko_error("aabb_contains: One or both aabbs are not initialized")
454  ! end if
455 
456  is_contained = all(this%box_min .le. other%box_min) .and. &
457  all(this%box_max .ge. other%box_max)
458 
459  end function aabb_contains_other
460 
462  function aabb_contains_point(this, p) result(is_contained)
463  class(aabb_t), intent(in) :: this
464  real(kind=dp), dimension(3), intent(in) :: p
465  logical :: is_contained
466 
467  ! if (.not. this%initialized) then
468  ! call neko_error("aabb_contains_point: One or both aabbs are not initialized")
469  ! end if
470 
471  is_contained = all(p .ge. this%box_min) .and. all(p .le. this%box_max)
472  end function aabb_contains_point
473 
474  ! ========================================================================== !
475  ! Binary operations
476  ! ========================================================================== !
477 
479  function merge_aabb(box1, box2) result(merged)
480  class(aabb_t), intent(in) :: box1
481  class(aabb_t), intent(in) :: box2
482  type(aabb_t) :: merged
483 
484  real(kind=dp), dimension(3) :: box_min, box_max
485 
486  box_min = min(box1%box_min, box2%box_min)
487  box_max = max(box1%box_max, box2%box_max)
488 
489  call merged%init(box_min, box_max)
490  end function merge_aabb
491 
493  function intersection_aabb(box1, box2) result(intersected)
494  class(aabb_t), intent(in) :: box1
495  class(aabb_t), intent(in) :: box2
496  type(aabb_t) :: intersected
497 
498  real(kind=dp), dimension(3) :: box_min, box_max
499 
500  box_min = max(box1%box_min, box2%box_min)
501  box_max = min(box1%box_max, box2%box_max)
502 
503  call intersected%init(box_min, box_max)
504  end function intersection_aabb
505 
506  ! ========================================================================== !
507  ! Private operations
508  ! ========================================================================== !
509 
511  pure function calculate_surface_area(this) result(surface_area)
512  class(aabb_t), intent(in) :: this
513  real(kind=dp) :: surface_area
514 
515  surface_area = 2.0 * (&
516  & this%get_width() * this%get_height() &
517  & + this%get_width() * this%get_depth() &
518  & + this%get_height() * this%get_depth() &
519  &)
520  end function calculate_surface_area
521 
522  ! ========================================================================== !
523  ! Comparison operators
524  ! ========================================================================== !
525 
527  pure function aabb_less(this, other)
528  class(aabb_t), intent(in) :: this
529  class(aabb_t), intent(in) :: other
530  logical :: aabb_less
531  logical :: equal
532 
533  if (.not. this%initialized .or. .not. other%initialized) then
534  aabb_less = .false.
535  return
536  end if
537 
538  aabb_less = this%box_min(1) .lt. other%box_min(1)
539  equal = this%box_min(1) .le. other%box_min(1)
540 
541  if (.not. aabb_less .and. equal) then
542  aabb_less = this%box_min(2) .lt. other%box_min(2)
543  equal = this%box_min(2) .le. other%box_min(2)
544  end if
545  if (.not. aabb_less .and. equal) then
546  aabb_less = this%box_min(3) .lt. other%box_min(3)
547  equal = this%box_min(3) .le. other%box_min(3)
548  end if
549  if (.not. aabb_less .and. equal) then
550  aabb_less = this%box_max(1) .lt. other%box_max(1)
551  equal = this%box_max(1) .le. other%box_max(1)
552  end if
553  if (.not. aabb_less .and. equal) then
554  aabb_less = this%box_max(2) .lt. other%box_max(2)
555  equal = this%box_max(2) .le. other%box_max(2)
556  end if
557  if (.not. aabb_less .and. equal) then
558  aabb_less = this%box_max(3) .lt. other%box_max(3)
559  end if
560 
561  end function aabb_less
562 
564  pure function aabb_greater(this, other)
565  class(aabb_t), intent(in) :: this
566  class(aabb_t), intent(in) :: other
567  logical :: aabb_greater
568  logical :: equal
569 
570  if (.not. this%initialized .or. .not. other%initialized) then
571  aabb_greater = .false.
572  return
573  end if
574 
575  aabb_greater = this%box_min(1) .gt. other%box_min(1)
576  equal = this%box_min(1) .ge. other%box_min(1)
577 
578  if (.not. aabb_greater .and. equal) then
579  aabb_greater = this%box_min(2) .gt. other%box_min(2)
580  equal = this%box_min(2) .ge. other%box_min(2)
581  end if
582  if (.not. aabb_greater .and. equal) then
583  aabb_greater = this%box_min(3) .gt. other%box_min(3)
584  equal = this%box_min(3) .ge. other%box_min(3)
585  end if
586  if (.not. aabb_greater .and. equal) then
587  aabb_greater = this%box_max(1) .gt. other%box_max(1)
588  equal = this%box_max(1) .ge. other%box_max(1)
589  end if
590  if (.not. aabb_greater .and. equal) then
591  aabb_greater = this%box_max(2) .gt. other%box_max(2)
592  equal = this%box_max(2) .ge. other%box_max(2)
593  end if
594  if (.not. aabb_greater .and. equal) then
595  aabb_greater = this%box_max(3) .gt. other%box_max(3)
596  end if
597 
598  end function aabb_greater
599 
600 end module aabb
Intersect two aabbs.
Definition: aabb.f90:98
Merge two aabbs.
Definition: aabb.f90:93
Axis Aligned Bounding Box (aabb) implementation in Fortran.
Definition: aabb.f90:71
pure real(kind=dp) function, dimension(3) aabb_get_max(this)
Get the maximum point of the aabb.
Definition: aabb.f90:363
pure real(kind=dp) function aabb_get_surface_area(this)
Get the surface area of the aabb.
Definition: aabb.f90:403
pure real(kind=dp) function aabb_get_diameter(this)
Get the diameter length of the aabb.
Definition: aabb.f90:395
subroutine add_padding(this, padding)
Add padding to the aabb.
Definition: aabb.f90:222
logical function aabb_contains_point(this, p)
Check if this aabb contains a point.
Definition: aabb.f90:463
subroutine aabb_init(this, lower_left_front, upper_right_back)
Initialize the aabb.
Definition: aabb.f90:335
type(aabb_t) function get_aabb_tri_mesh(object)
Get the aabb of a triangular mesh.
Definition: aabb.f90:291
logical function aabb_overlaps(this, other)
Check if two aabbs are overlapping.
Definition: aabb.f90:431
type(aabb_t) function get_aabb_mesh(object)
Get the aabb of a mesh.
Definition: aabb.f90:269
type(aabb_t) function intersection_aabb(box1, box2)
Get the intersection of two aabbs.
Definition: aabb.f90:494
pure real(kind=dp) function aabb_get_height(this)
Get the height of the aabb. Also known as the z-axis length.
Definition: aabb.f90:387
pure logical function aabb_greater(this, other)
Greater than comparison operator.
Definition: aabb.f90:565
type(aabb_t) function, public get_aabb(object, padding)
Construct the aabb of a predefined object.
Definition: aabb.f90:181
logical function aabb_contains_other(this, other)
Check if this aabb contains another aabb.
Definition: aabb.f90:448
pure real(kind=dp) function calculate_surface_area(this)
Calculate the surface area of the aabb.
Definition: aabb.f90:512
pure logical function aabb_less(this, other)
Less than comparison operator.
Definition: aabb.f90:528
type(aabb_t) function get_aabb_element(object)
Get the aabb of an arbitrary element.
Definition: aabb.f90:242
type(aabb_t) function get_aabb_tet_mesh(object)
Get the aabb of a tetrahedral mesh.
Definition: aabb.f90:313
pure real(kind=dp) function, dimension(3) aabb_get_center(this)
Get the center of the aabb.
Definition: aabb.f90:411
pure real(kind=dp) function aabb_get_width(this)
Get the width of the aabb. Also known as the x-axis length.
Definition: aabb.f90:371
pure real(kind=dp) function, dimension(3) aabb_get_diagonal(this)
Get the diagonal of the aabb.
Definition: aabb.f90:419
type(aabb_t) function merge_aabb(box1, box2)
Merge two aabbs.
Definition: aabb.f90:480
pure real(kind=dp) function aabb_get_depth(this)
Get the depth of the aabb. Also known as the y-axis length.
Definition: aabb.f90:379
pure real(kind=dp) function, dimension(3) aabb_get_min(this)
Get the minimum point of the aabb.
Definition: aabb.f90:355
Defines a hexahedron element.
Definition: hex.f90:34
Defines a mesh.
Definition: mesh.f90:34
integer, parameter, public dp
Definition: num_types.f90:9
Implements a point.
Definition: point.f90:35
Defines a quadrilateral element.
Definition: quad.f90:34
Defines a tetrahedral mesh.
Definition: tet_mesh.f90:35
Defines a tetrahedral element.
Definition: tet.f90:34
Defines a triangular surface mesh.
Definition: tri_mesh.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
Base type for an element.
Definition: element.f90:44
Hexahedron element.
Definition: hex.f90:63
A point in with coordinates .
Definition: point.f90:43
Quadrilateral element.
Definition: quad.f90:58
Tetrahedral element.
Definition: tet.f90:63
Triangular element.
Definition: tri.f90:60
#define max(a, b)
Definition: tensor.cu:40