Neko  0.9.0
A portable framework for high-order spectral element flow simulations
octree.f90
Go to the documentation of this file.
1 ! Copyright (c) 2022, 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 !
35 module octree
36  use num_types
37  use point
38  use utils
39  implicit none
40  private
41 
42  type oct_ptr_t
43  type(oct_t), pointer :: ptr => null()
44  end type oct_ptr_t
45 
47  type oct_t
48  type(oct_ptr_t) :: oct(8)
49  integer :: level
50  type(point_t) :: origin
51  real(kind=dp) :: width
52  type(point_t) :: point
53  logical valid
54  end type oct_t
55 
57  type, public :: octree_t
58  type(oct_t), pointer :: root => null()
59  contains
60  procedure, pass(t) :: init => octree_init
61  procedure, pass(t) :: free => octree_free
62  procedure, pass(t) :: insert => octree_insert
63  procedure, pass(t) :: find => octree_find
64  end type octree_t
65 
66 contains
67 
69  subroutine octree_init(t, width)
70  class(octree_t), intent(inout) :: t
71  real(kind=dp), intent(in) :: width
72  type(point_t) :: origin
73  integer, parameter :: top_level = 0
74 
75  call octree_free(t)
76 
77  origin = (/ 0d0, 0d0, 0d0 /)
78  call octree_oct_init(t%root, origin, width, top_level)
79 
80  end subroutine octree_init
81 
83  subroutine octree_free(t)
84  class(octree_t), intent(inout) :: t
85 
86  call octree_free_oct(t%root)
87 
88  end subroutine octree_free
89 
91  subroutine octree_insert(t, p)
92  class(octree_t), intent(inout) :: t
93  type(point_t), intent(in) :: p
94 
95  call octree_oct_insert(t%root, p)
96 
97  end subroutine octree_insert
98 
100  function octree_find(t, p) result(rcode)
101  class(octree_t), intent(in) :: t
102  type(point_t), intent(in) :: p
103  logical rcode
104 
105  rcode = (octree_oct_find(t%root, p) .eq. 0)
106 
107  end function octree_find
108 
109 
111  recursive subroutine octree_oct_insert(o, p)
112  type(oct_t), intent(inout), pointer :: o
113  type(point_t), intent(in) :: p
114  type(point_t) :: tmp_pt, offset, new_origin
115  integer :: i
116 
117  if (.not. associated(o%oct(1)%ptr)) then
118  if (.not. o%valid) then
119  o%point = p
120  o%valid = .true.
121  return
122  else
123  if (o%point .eq. p) then
124  return
125  else
126  tmp_pt = o%point
127  o%valid = .false.
128 
129  do i = 1, 8
130  offset = (/ -0.5d0, -0.5d0, -0.5d0 /)
131  if (iand((i - 1), 4) .gt. 0) offset%x(1) = 0.5d0
132  if (iand((i - 1), 2) .gt. 0) offset%x(2) = 0.5d0
133  if (iand((i - 1), 1) .gt. 0) offset%x(3) = 0.5d0
134 
135  new_origin = o%origin%x + (o%width * offset%x)
136  call octree_oct_init(o%oct(i)%ptr, new_origin, &
137  o%width * 0.5d0, o%level + 1)
138  end do
139 
140  call octree_oct_insert(o%oct(octree_oct(o, tmp_pt))%ptr, tmp_pt)
141  call octree_oct_insert(o%oct(octree_oct(o, p))%ptr, p)
142 
143  end if
144  end if
145  else
146  call octree_oct_insert(o%oct(octree_oct(o, p))%ptr, p)
147  end if
148  end subroutine octree_oct_insert
149 
151  recursive function octree_oct_find(o, p) result(rcode)
152  type(oct_t), pointer, intent(in) :: o
153  type(point_t), intent(in) :: p
154  integer :: rcode
155  integer :: oct_idx
156 
157  rcode = 1
158 
159  if (.not. associated(o%oct(1)%ptr)) then
160  if (o%valid .and. octree_oct_inside(o, p)) then
161  rcode = 0
162  return
163  end if
164  else
165  oct_idx = octree_oct(o, p)
166  rcode = octree_oct_find(o%oct(oct_idx)%ptr, p)
167  end if
168 
169  end function octree_oct_find
170 
172  subroutine octree_oct_init(o, origin, width, level)
173  type(oct_t), pointer, intent(inout) :: o
174  type(point_t), intent(in) :: origin
175  real(kind=dp), intent(in) :: width
176  integer, intent(in) :: level
177  integer :: i
178 
179  if (associated(o)) then
180  call neko_error('Octree octant already initialized')
181  else
182  allocate(o)
183  o%origin = origin
184  o%width = width
185  o%level = level
186  o%valid = .false.
187 
188  do i = 1, 8
189  nullify(o%oct(i)%ptr)
190  end do
191  end if
192 
193  end subroutine octree_oct_init
194 
196  recursive subroutine octree_free_oct(o)
197  type(oct_t), pointer, intent(inout) :: o
198  integer :: i
199 
200  if (.not. associated(o)) then
201  return
202  else if (.not. associated(o%oct(1)%ptr)) then
203  deallocate(o)
204  nullify(o)
205  return
206  else
207  do i = 1, 8
208  call octree_free_oct(o%oct(i)%ptr)
209  end do
210  deallocate(o)
211  nullify(o)
212  end if
213 
214 
215  end subroutine octree_free_oct
216 
218  pure function octree_oct(oct, point) result(oct_idx)
219  type(oct_t), pointer, intent(in) :: oct
220  type(point_t), intent(in) :: point
221  integer :: oct_idx
222 
223  oct_idx = 0
224  if (point%x(1) .ge. oct%origin%x(1)) oct_idx = ior(oct_idx, 4)
225  if (point%x(2) .ge. oct%origin%x(2)) oct_idx = ior(oct_idx, 2)
226  if (point%x(3) .ge. oct%origin%x(3)) oct_idx = ior(oct_idx, 1)
227  oct_idx = oct_idx + 1
228 
229  end function octree_oct
230 
232  pure function octree_oct_inside(oct, point) result(inside)
233  type(oct_t), pointer, intent(in) :: oct
234  type(point_t), intent(in) :: point
235  logical :: inside
236 
237  inside = .false.
238  if ((point%x(1) .ge. (oct%origin%x(1) - oct%width) .and. &
239  point%x(1) .le. (oct%origin%x(1) + oct%width)) .and. &
240  (point%x(2) .ge. (oct%origin%x(2) - oct%width) .and. &
241  point%x(2) .le. (oct%origin%x(2) + oct%width)) .and. &
242  (point%x(3) .ge. (oct%origin%x(3) - oct%width) .and. &
243  point%x(3) .le. (oct%origin%x(3) + oct%width))) then
244  inside = .true.
245  end if
246  end function octree_oct_inside
247 
248 end module octree
integer, parameter, public dp
Definition: num_types.f90:9
Implements an Octree.
Definition: octree.f90:35
pure logical function octree_oct_inside(oct, point)
Return if a point is inside an octant.
Definition: octree.f90:233
subroutine octree_free(t)
Destroy an octree.
Definition: octree.f90:84
subroutine octree_insert(t, p)
Insert a point p into the octree.
Definition: octree.f90:92
recursive subroutine octree_oct_insert(o, p)
Insert a point p into the octree rooted at o.
Definition: octree.f90:112
recursive subroutine octree_free_oct(o)
Deallocate an oct in an octree.
Definition: octree.f90:197
logical function octree_find(t, p)
Find a point p in an octree.
Definition: octree.f90:101
pure integer function octree_oct(oct, point)
Return the octant for a given point.
Definition: octree.f90:219
subroutine octree_oct_init(o, origin, width, level)
Initialize an octant width a given width, origin and level.
Definition: octree.f90:173
recursive integer function octree_oct_find(o, p)
Find the octant containing a point p.
Definition: octree.f90:152
subroutine octree_init(t, width)
Initialize an octree.
Definition: octree.f90:70
Implements a point.
Definition: point.f90:35
Utilities.
Definition: utils.f90:35
Defines an octree octant.
Definition: octree.f90:47
Defines an octree.
Definition: octree.f90:57
A point in with coordinates .
Definition: point.f90:43