File: N:\mfix\model\des\sap.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module: sweep_and_prune                                             !
4     !                                                                      !
5     !  Purpose: maintains a sorted list of axis-aligned bounding boxes     !
6     !                                                                      !
7     !           main type is sap_t                                         !
8     !                                                                      !
9     !  Reference: http://www.codercorner.com/SAP.pdf                       !
10     !                                                                      !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
12     
13     #include "version.inc"
14     
15     module sweep_and_prune
16     
17       use pair_manager
18       use list
19     
20       ! endpoint data structure should be 64 bits in size
21       type endpoint_t
22          ! owner of the endpoint
23          ! if Min endpoint, id = -box_id
24          ! if Max endpoint, id = box_id
25          integer(kind=4) :: box_id
26     
27          ! actual value of the endpoint; single precision to save space
28          real :: value
29       end type endpoint_t
30     
31       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
32       !                                                                      !
33       !  Type: aabb_t                                                        !
34       !                                                                      !
35       !  Purpose: Represents an axis-aligned bounding box                    !
36       !           Used to add new boxes to sap_t and multisap_t              !
37       !                                                                      !
38       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
39     
40       type aabb_t
41          real :: maxendpoint(3)
42          real :: minendpoint(3)
43       end type aabb_t
44     
45       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
46       !                                                                      !
47       !  Type: box_t                                                         !
48       !                                                                      !
49       !  Purpose: Represents the AABB of a particle                          !
50       !           maxendpoint(1) and minendpoint_id(1) point to x_endpoints  !
51       !           maxendpoint(2) and minendpoint_id(2) point to y_endpoints  !
52       !           maxendpoint(3) and minendpoint_id(3) point to z_endpoints  !
53       !                                                                      !
54       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
55     
56       type box_t
57          integer :: maxendpoint_id(3)
58          integer :: minendpoint_id(3)
59          integer :: particle_id
60       end type box_t
61     
62       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
63       !                                                                      !
64       !  Type: sap_t                                                         !
65       !                                                                      !
66       !  Purpose: Sorts the endpoints array                                  !
67       !           axis is (1,2,3) to represent x,y,z                         !
68       !                                                                      !
69       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
70     
71       type sap_t
72          ! id to identify within a multisap
73          integer :: id
74     
75          ! list of endpoints of boxes
76          type(endpoint_t), dimension(:), allocatable :: x_endpoints
77          type(endpoint_t), dimension(:), allocatable :: y_endpoints
78          type(endpoint_t), dimension(:), allocatable :: z_endpoints
79          integer :: x_endpoints_len
80          integer :: y_endpoints_len
81          integer :: z_endpoints_len
82     
83          ! list of boxes
84          type(box_t), dimension(:), allocatable :: boxes
85          integer :: boxes_len
86     
87          ! list of indices in boxes that have value 0 (deleted boxes)
88          type(list_t) :: deleted_boxes
89     
90          ! list of pairs of particle ids that have intersecting boxes
91          type(hashtable_t) :: hashtable
92       end type sap_t
93     
94       public :: init_sap, add_box, del_box, update_box, sort, sweep
95       private :: partition, fullcheck, boxes_grow, endpoints_grow, sort_endpoints, quicksort_endpoints
96     
97     contains
98     
99       subroutine check(this)
100         use discretelement
101         implicit none
102         type(sap_t), intent(inout) :: this
103         integer :: ii, pid
104         do ii = 1, this%x_endpoints_len
105     
106            if (0 > this%x_endpoints(ii)%box_id) then
107     
108            if (ii .ne. this%boxes(-this%x_endpoints(ii)%box_id)%minendpoint_id(1)  ) then
109               print *,"ii is ",ii
110               print *,"this%x_endpoints(ii)%box_id is ",this%x_endpoints(ii)%box_id
111               print *,"this%boxes(this%x_endpoints(ii)%box_id)%minendpoint_id(1)",this%boxes(-this%x_endpoints(ii)%box_id)%minendpoint_id(1)
112               stop __LINE__
113            endif
114     else
115            if (ii .ne. this%boxes(this%x_endpoints(ii)%box_id)%maxendpoint_id(1)  ) then
116               print *,"ii is ",ii
117               print *,"this%x_endpoints(ii)%box_id is ",this%x_endpoints(ii)%box_id
118               print *,"this%boxes(this%x_endpoints(ii)%box_id)%minendpoint_id(1)",this%boxes(this%x_endpoints(ii)%box_id)%maxendpoint_id(1)
119     
120               stop __LINE__
121            endif
122     endif
123         enddo
124     
125         do ii = 1, this%y_endpoints_len
126            if (0 > this%y_endpoints(ii)%box_id) then
127               if (ii .ne. this%boxes(-this%y_endpoints(ii)%box_id)%minendpoint_id(2)  ) then
128                  print *,"ii is ",ii
129                  print *,"this%y_endpoints(ii)%box_id is ",this%y_endpoints(ii)%box_id
130                  print *,"this%boxes(this%y_endpoints(ii)%box_id)%minendpoint_id(1)",this%boxes(-this%y_endpoints(ii)%box_id)%minendpoint_id(2)
131     
132                  stop __LINE__
133               endif
134            else
135               if (ii .ne. this%boxes(this%y_endpoints(ii)%box_id)%maxendpoint_id(2)  ) then
136                  stop __LINE__
137               endif
138            endif
139         enddo
140     
141         do ii = 1, this%z_endpoints_len
142            if (0 > this%z_endpoints(ii)%box_id) then
143               if (ii .ne. this%boxes(-this%z_endpoints(ii)%box_id)%minendpoint_id(3)  ) then
144                  stop __LINE__
145               endif
146            else
147               if (ii .ne. this%boxes(this%z_endpoints(ii)%box_id)%maxendpoint_id(3)  ) then
148                  stop __LINE__
149               endif
150            endif
151         enddo
152     
153         do ii = 1, this%boxes_len
154            pid = this%boxes(ii)%particle_id
155            if ( 0.0001 < abs(this%x_endpoints(this%boxes(ii)%maxendpoint_id(1))%value - this%x_endpoints(this%boxes(ii)%minendpoint_id(1))%value - 2*des_radius(pid))) then
156     
157               print *,this%x_endpoints(this%boxes(ii)%maxendpoint_id(1))%value
158               print *,this%x_endpoints(this%boxes(ii)%minendpoint_id(1))%value
159               print *,2*des_radius(pid)
160     
161               stop __LINE__
162            endif
163         enddo
164     
165         do ii = 1, this%boxes_len
166            pid = this%boxes(ii)%particle_id
167            if ( 0.0001 < abs(this%y_endpoints(this%boxes(ii)%maxendpoint_id(2))%value - this%y_endpoints(this%boxes(ii)%minendpoint_id(2))%value - 2*des_radius(pid))) then
168               stop __LINE__
169            endif
170         enddo
171     
172         do ii = 1, this%boxes_len
173            pid = this%boxes(ii)%particle_id
174            if ( 0.0001 < abs(this%z_endpoints(this%boxes(ii)%maxendpoint_id(3))%value - this%z_endpoints(this%boxes(ii)%minendpoint_id(3))%value - 2*des_radius(pid))) then
175               stop __LINE__
176            endif
177         enddo
178     
179       end subroutine check
180     
181       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
182       !                                                                      !
183       !  Subroutine: init_sap                                                !
184       !                                                                      !
185       !  Purpose: sap_t constructor                                          !
186       !                                                                      !
187       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
188     
189       subroutine init_sap(this,id)
190         implicit none
191         type(sap_t), intent(inout) :: this
192         integer, intent(in) :: id
193     
194         this%x_endpoints_len = 0
195         this%y_endpoints_len = 0
196         this%z_endpoints_len = 0
197         this%boxes_len = 0
198         this%id = id
199     
200         allocate (this%x_endpoints(10))
201         allocate (this%y_endpoints(10))
202         allocate (this%z_endpoints(10))
203         allocate (this%boxes(10))
204     
205         call list_init(this%deleted_boxes)
206     
207         call init_pairs(this%hashtable)
208     
209       end subroutine init_sap
210     
211       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
212       !                                                                      !
213       !  Subroutine: boxes_GROW                                              !
214       !                                                                      !
215       !  Purpose: resize array of box_t                                      !
216       !                                                                      !
217       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
218     
219       SUBROUTINE boxes_GROW(box_array,new_size)
220         IMPLICIT NONE
221     
222         INTEGER, INTENT(IN) :: new_size
223         type(box_t), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: box_array
224         type(box_t), DIMENSION(:), ALLOCATABLE :: box_tmp
225         INTEGER lSIZE
226     
227         lSIZE = size(box_array,1)
228         allocate(box_tmp(new_size))
229         box_tmp(1:lSIZE) = box_array(1:lSIZE)
230         call move_alloc(box_tmp,box_array)
231     
232       END SUBROUTINE boxes_GROW
233     
234       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
235       !                                                                      !
236       !  Subroutine: endpoints_GROW                                          !
237       !                                                                      !
238       !  Purpose: Resizes array of endpoint_t                                !
239       !                                                                      !
240       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
241     
242       SUBROUTINE endpoints_GROW(endpoint_array,new_size)
243         IMPLICIT NONE
244     
245         INTEGER, INTENT(IN) :: new_size
246         type(endpoint_t), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: endpoint_array
247         type(endpoint_t), DIMENSION(:), ALLOCATABLE :: endpoint_tmp
248         INTEGER lSIZE
249     
250         lSIZE = size(endpoint_array,1)
251         allocate(endpoint_tmp(new_size))
252         endpoint_tmp(1:lSIZE) = endpoint_array(1:lSIZE)
253         call move_alloc(endpoint_tmp,endpoint_array)
254     
255       END SUBROUTINE endpoints_GROW
256     
257       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
258       !                                                                      !
259       !  Subroutine: add_box                                                 !
260       !                                                                      !
261       !  Purpose: Adds new box_t to sap_t that represents particle_id        !
262       !           Returns id to refer to the box later                       !
263       !                                                                      !
264       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
265     
266       subroutine add_box(this,aabb,particle_id,id)
267     
268         implicit none
269         type(sap_t), intent(inout) :: this
270         type(aabb_t), intent(in) :: aabb
271         integer, intent(in) :: particle_id
272         integer, intent(out) :: id
273     
274         if (this%x_endpoints_len+2 > size(this%x_endpoints)) then
275            call endpoints_GROW(this%x_endpoints,2*(this%x_endpoints_len+2))
276         endif
277     
278         if (this%y_endpoints_len+2 > size(this%y_endpoints)) then
279            call endpoints_GROW(this%y_endpoints,2*(this%y_endpoints_len+2))
280         endif
281     
282         if (this%z_endpoints_len+2 > size(this%z_endpoints)) then
283            call endpoints_GROW(this%z_endpoints,2*(this%z_endpoints_len+2))
284         endif
285     
286         if (this%boxes_len+1 > size(this%boxes)) then
287            call boxes_GROW(this%boxes,2*(this%boxes_len+1))
288         endif
289     
290         id = list_pop(this%deleted_boxes)
291     
292         ! if there are no deleted boxes (no holes in boxes array), then append the
293         ! new box at the end of the boxes array.
294         if ( id < 0 ) then
295            this%boxes_len = this%boxes_len + 1
296            id = this%boxes_len
297         endif
298     
299         ! no holes in the endpoints array; new endpoints are always appended
300         this%x_endpoints_len = this%x_endpoints_len + 2
301         this%y_endpoints_len = this%y_endpoints_len + 2
302         this%z_endpoints_len = this%z_endpoints_len + 2
303     
304         this%boxes(id)%particle_id = particle_id
305         this%boxes(id)%maxendpoint_id(1) = this%x_endpoints_len
306         this%boxes(id)%maxendpoint_id(2) = this%y_endpoints_len
307         this%boxes(id)%maxendpoint_id(3) = this%z_endpoints_len
308         this%boxes(id)%minendpoint_id(1) = this%x_endpoints_len-1
309         this%boxes(id)%minendpoint_id(2) = this%y_endpoints_len-1
310         this%boxes(id)%minendpoint_id(3) = this%z_endpoints_len-1
311     
312         this%x_endpoints(this%x_endpoints_len-1)%box_id = -id
313         this%x_endpoints(this%x_endpoints_len-1)%value = aabb%minendpoint(1)
314         this%y_endpoints(this%y_endpoints_len-1)%box_id = -id
315         this%y_endpoints(this%y_endpoints_len-1)%value = aabb%minendpoint(2)
316         this%z_endpoints(this%z_endpoints_len-1)%box_id = -id
317         this%z_endpoints(this%z_endpoints_len-1)%value = aabb%minendpoint(3)
318     
319         this%x_endpoints(this%x_endpoints_len)%box_id = id
320         this%x_endpoints(this%x_endpoints_len)%value = aabb%maxendpoint(1)
321         this%y_endpoints(this%y_endpoints_len)%box_id = id
322         this%y_endpoints(this%y_endpoints_len)%value = aabb%maxendpoint(2)
323         this%z_endpoints(this%z_endpoints_len)%box_id = id
324         this%z_endpoints(this%z_endpoints_len)%value = aabb%maxendpoint(3)
325     
326       end subroutine add_box
327     
328       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
329       !                                                                      !
330       !  Subroutine: del_box                                                 !
331       !                                                                      !
332       !  Purpose: Removes box from sap_t                                     !
333       !           Setting all the endpoints to HUGE will remove all pairs    !
334       !           associated with this box as the endpoints get sorted.      !
335       !           The box is "really" deleted at the end of sort().          !
336       !                                                                      !
337       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
338     
339       subroutine del_box(this,id)
340         implicit none
341         type(sap_t), intent(inout) :: this
342         integer, intent(in) :: id
343     
344         this%x_endpoints(this%boxes(id)%maxendpoint_id(1))%value = HUGE(0.0)
345         this%y_endpoints(this%boxes(id)%maxendpoint_id(2))%value = HUGE(0.0)
346         this%z_endpoints(this%boxes(id)%maxendpoint_id(3))%value = HUGE(0.0)
347         this%x_endpoints(this%boxes(id)%minendpoint_id(1))%value = HUGE(0.0)
348         this%y_endpoints(this%boxes(id)%minendpoint_id(2))%value = HUGE(0.0)
349         this%z_endpoints(this%boxes(id)%minendpoint_id(3))%value = HUGE(0.0)
350     
351       end subroutine del_box
352     
353       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
354       !                                                                      !
355       !  Subroutine: update_box                                              !
356       !                                                                      !
357       !  Purpose: Update box represented by id with values in aabb           !
358       !                                                                      !
359       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
360     
361       subroutine update_box(this,id,aabb)
362         implicit none
363         type(sap_t), intent(inout) :: this
364         integer, intent(in) :: id
365         type(aabb_t), intent(in) :: aabb
366     
367         ! update end points
368         this%x_endpoints(this%boxes(id)%maxendpoint_id(1))%value = aabb%maxendpoint(1)
369         this%y_endpoints(this%boxes(id)%maxendpoint_id(2))%value = aabb%maxendpoint(2)
370         this%z_endpoints(this%boxes(id)%maxendpoint_id(3))%value = aabb%maxendpoint(3)
371         this%x_endpoints(this%boxes(id)%minendpoint_id(1))%value = aabb%minendpoint(1)
372         this%y_endpoints(this%boxes(id)%minendpoint_id(2))%value = aabb%minendpoint(2)
373         this%z_endpoints(this%boxes(id)%minendpoint_id(3))%value = aabb%minendpoint(3)
374     
375       end subroutine update_box
376     
377       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
378       !                                                                      !
379       !  Subroutine: sort                                                    !
380       !                                                                      !
381       !  Purpose: (Re)sorts particles after each timestep,                   !
382       !           which keeps the pair manager updated.                      !
383       !           Calls sort_endpoints, which does an insertion sort,        !
384       !           which is fast for nearly sorted arrays but slow for        !
385       !           unsorted arrays.                                           !
386       !                                                                      !
387       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
388     
389       subroutine sort(this)
390         use geometry, only: do_k
391         implicit none
392         type(sap_t), intent(inout) :: this
393         integer :: ii
394     
395         ! sort end points
396         call sort_endpoints(this,this%x_endpoints(1:this%x_endpoints_len),1)
397         call sort_endpoints(this,this%y_endpoints(1:this%y_endpoints_len),2)
398         if (do_k) call sort_endpoints(this,this%z_endpoints(1:this%z_endpoints_len),3)
399     
400         ! cleanup HUGE endpoints
401         !
402         ! After sorting the endpoints, there may be endpoints at the end of the
403         ! array corresponding to boxes that are being deleted. Now that the sorting
404         ! is done, the corresponding ids in boxes array can be set to zero, and
405         ! added to the deleted boxes array to be used next time we add a box.
406     
407         ! TODO: what if two overlapping boxes are deleted at the same time...could
408         ! they leave a pair in the pair manager corresponding to deleted particles?
409         ! Not sure how to handle that
410     
411         ii = this%x_endpoints_len
412         do
413            if (ii .le. 0) exit
414            if (this%x_endpoints(ii)%value .ne. HUGE(0.0)) exit
415            this%boxes(abs(this%x_endpoints(ii)%box_id))%minendpoint_id(:) = 0
416            this%boxes(abs(this%x_endpoints(ii)%box_id))%maxendpoint_id(:) = 0
417            call list_add(this%deleted_boxes,this%x_endpoints(ii)%box_id)
418            ii = ii - 1
419         enddo
420     
421         this%x_endpoints_len = ii
422         this%y_endpoints_len = ii
423         this%z_endpoints_len = ii
424     
425       end subroutine sort
426     
427       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
428       !                                                                      !
429       !  Subroutine: sweep                                                   !
430       !                                                                      !
431       !  Purpose: Creates initial set of pairs of colliding boxes            !
432       !           Called on startup after quicksort is called                !
433       !           for each axis set of endpoints                             !
434       !                                                                      !
435       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
436     
437       subroutine sweep(this)
438     
439         use pair_manager, only: add_pair
440         use geometry, only: no_k
441     
442         implicit none
443         type(sap_t), intent(inout) :: this
444     
445         ! active list of box id's
446         integer :: ii,aa,minmax,ai
447     
448         type(list_t) :: active
449     
450         call list_init(active)
451     
452         do ii=1, this%x_endpoints_len
453     
454            minmax = this%x_endpoints(ii)%box_id
455     
456            if ( minmax < 0) then
457               ! add pairs for new box x ( active pairs )
458     
459               do ai=1, list_get_length(active)
460                  aa = list_get(active,ai)
461                  ! compare other two axes
462     
463                  if (max(this%boxes(aa)%minendpoint_id(2),this%boxes(-minmax)%minendpoint_id(2)) <= min(this%boxes(-minmax)%maxendpoint_id(2),this%boxes(aa)%maxendpoint_id(2)) .and. &
464                       (NO_K .or. max(this%boxes(aa)%minendpoint_id(3),this%boxes(-minmax)%minendpoint_id(3)) <= min(this%boxes(-minmax)%maxendpoint_id(3),this%boxes(aa)%maxendpoint_id(3)))) then
465                     call add_pair(this%hashtable,this%boxes(-minmax)%particle_id,this%boxes(aa)%particle_id)
466                  endif
467               enddo
468     
469               ! add new box to active pair list
470               call list_add(active,-minmax)
471     
472            else if ( 0 < minmax ) then
473               ! remove box from active pair list
474               call list_del(active,minmax)
475            else
476               print *,"SAP_ID=",this%id,"minmax shouldn't be zero: ",minmax
477               ERROR_STOP __LINE__
478            endif
479         enddo
480     
481       end subroutine sweep
482     
483       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
484       !                                                                      !
485       !  Subroutine: sort_endpoints                                          !
486       !                                                                      !
487       !  Purpose: Sorts the endpoints array                                  !
488       !           axis is (1,2,3) to represent x,y,z                         !
489       !                                                                      !
490       !           Sorts endpoints array after each timestep.                 !
491       !           Uses insertion sort.                                       !
492       !                                                                      !
493       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
494     
495       subroutine sort_endpoints(this, endpoints, axis)
496         use pair_manager, only: add_pair, del_pair
497         implicit none
498         type(sap_t), intent(inout) :: this
499         integer, intent(in) :: axis
500         type(endpoint_t), dimension(:), intent(inout) :: endpoints
501         integer :: ii, jj
502         type(endpoint_t) :: sweeppoint, swappoint
503         real :: sweepval
504     
505         do ii=2, size(endpoints)
506     
507            sweeppoint = endpoints(ii)
508            sweepval = sweeppoint%value
509            jj = ii-1
510     
511            do while ( 0 < jj )
512               if ( endpoints(jj)%value <= sweepval ) then
513                  exit
514               endif
515     
516               swappoint = endpoints(jj)
517     
518               if (sweeppoint%box_id < 0 .and. 0 < swappoint%box_id ) then
519                  if (fullcheck(this,-sweeppoint%box_id,swappoint%box_id,axis)) then
520                     call add_pair(this%hashtable,this%boxes(-sweeppoint%box_id)%particle_id,this%boxes(swappoint%box_id)%particle_id)
521                  endif
522               endif
523     
524               if (0 < sweeppoint%box_id .and. swappoint%box_id < 0 ) then
525                  if (fullcheck(this,sweeppoint%box_id,-swappoint%box_id,axis)) then
526                     call del_pair(this%hashtable,this%boxes(sweeppoint%box_id)%particle_id,this%boxes(-swappoint%box_id)%particle_id)
527                  endif
528               endif
529     
530               endpoints(jj+1) = swappoint
531     
532               ! if box_id is zero, that is an endpoint corresponding to a deleted box
533               if (swappoint%box_id < 0) then
534                  this%boxes(-swappoint%box_id)%minendpoint_id(axis) = jj+1
535               else if (swappoint%box_id > 0) then
536                  this%boxes(swappoint%box_id)%maxendpoint_id(axis) = jj+1
537               endif
538     
539               jj = jj - 1
540            enddo
541     
542            endpoints(jj+1) = sweeppoint
543     
544            ! if box_id is zero, that is an endpoint corresponding to a deleted box
545            if (sweeppoint%box_id < 0) then
546               this%boxes(-sweeppoint%box_id)%minendpoint_id(axis) = jj+1
547            else if (sweeppoint%box_id > 0) then
548               this%boxes(sweeppoint%box_id)%maxendpoint_id(axis) = jj+1
549            endif
550     
551         enddo
552     
553       end subroutine sort_endpoints
554     
555       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
556       !                                                                      !
557       !  Function: fullcheck                                                 !
558       !                                                                      !
559       !  Purpose: Check whether or checking that the boxes for id and id2    !
560       !           are correctly sorted                                       !
561       !                                                                      !
562       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
563     
564       logical function fullcheck(this,id,id2,curr_axis)
565         use geometry, only: no_k
566         implicit none
567         type(sap_t), intent(inout) :: this
568         ! box ids to compare
569         integer, intent(in) :: id, id2, curr_axis
570     
571         fullcheck = ((curr_axis.eq.1 .or. max(this%boxes(id)%minendpoint_id(1),this%boxes(id2)%minendpoint_id(1)) &
572                                           <= min(this%boxes(id2)%maxendpoint_id(1),this%boxes(id)%maxendpoint_id(1))) &
573                .and. (curr_axis.eq.2 .or. max(this%boxes(id)%minendpoint_id(2),this%boxes(id2)%minendpoint_id(2)) &
574                                            <= min(this%boxes(id2)%maxendpoint_id(2),this%boxes(id)%maxendpoint_id(2))) &
575      .and. (NO_K .or. curr_axis.eq.3 .or. max(this%boxes(id)%minendpoint_id(3),this%boxes(id2)%minendpoint_id(3)) &
576                                            <= min(this%boxes(id2)%maxendpoint_id(3),this%boxes(id)%maxendpoint_id(3))))
577     
578       end function fullcheck
579     
580       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
581       !                                                                      !
582       !  Subroutine: quicksort_endpoints                                     !
583       !                                                                      !
584       !  Purpose: Sorts the endpoints arrays                                 !
585       !                                                                      !
586       !           Only used once at startup, just before a call to sweep.    !
587       !           Quicksort is fast for unsorted arrays, but slow for        !
588       !           nearly sorted arrays.  sort() is used instead after each   !
589       !           timestep.                                                  !
590       !                                                                      !
591       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
592     
593       subroutine quicksort(this)
594         implicit none
595         type(sap_t), intent(inout) :: this
596     
597         call quicksort_endpoints(this,this%x_endpoints(1:this%x_endpoints_len),1)
598         call quicksort_endpoints(this,this%y_endpoints(1:this%y_endpoints_len),2)
599         call quicksort_endpoints(this,this%z_endpoints(1:this%z_endpoints_len),3)
600     
601       end subroutine quicksort
602     
603       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
604       !                                                                      !
605       !  Subroutine: quicksort_endpoints                                     !
606       !                                                                      !
607       !  Purpose: Sorts the endpoints array                                  !
608       !           axis is (1,2,3) to represent x,y,z                         !
609       !                                                                      !
610       !           Only used once at startup, just before a call to sweep.    !
611       !           Quicksort is fast for unsorted arrays, but slow for        !
612       !           nearly sorted arrays.  sort() is used instead after each   !
613       !           timestep.                                                  !
614       !                                                                      !
615       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
616     
617       recursive subroutine quicksort_endpoints(this, endpoints, axis)
618         type(sap_t), intent(inout) :: this
619         integer, intent(in) :: axis
620         type(endpoint_t), intent(inout), dimension(:) :: endpoints
621         integer :: ii
622     
623         !call check_boxes(this)
624     
625         if(size(endpoints) > 1) then
626            call partition(this, endpoints, ii, axis)
627            call quicksort_endpoints(this,endpoints(:ii-1),axis)
628            call quicksort_endpoints(this,endpoints(ii:),axis)
629         endif
630     
631       end subroutine quicksort_endpoints
632     
633       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
634       !                                                                      !
635       !  subroutine: partition                                               !
636       !                                                                      !
637       !  purpose: partition subroutine for quicksort of endpoints array      !
638       !           box is represented by id                                   !
639       !                                                                      !
640       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
641     
642       subroutine partition(this, endpoints, marker, axis)
643         type(sap_t), intent(inout) :: this
644         integer, intent(in) :: axis
645         type(endpoint_t), intent(in out), dimension(:) :: endpoints
646         integer, intent(out) :: marker
647         type(endpoint_t) :: temp
648         integer :: tmp_ii, tmp_jj
649         integer :: i, j
650         real :: x      ! pivot point
651         x = endpoints(1)%value
652         i= 0
653         j= size(endpoints) + 1
654     
655         do
656            j = j-1
657            do
658               if (endpoints(j)%value <= x) exit
659               j = j-1
660            end do
661            i = i+1
662            do
663               if (endpoints(i)%value >= x) exit
664               i = i+1
665            end do
666            if (i < j) then
667               ! exchange endpoints(i) and endpoints(j)
668     
669               if (endpoints(i)%box_id < 0) then
670                  tmp_ii = this%boxes(-endpoints(i)%box_id)%minendpoint_id(axis)
671               else
672                  tmp_ii = this%boxes(endpoints(i)%box_id)%maxendpoint_id(axis)
673               endif
674               if (endpoints(j)%box_id < 0) then
675                  tmp_jj = this%boxes(-endpoints(j)%box_id)%minendpoint_id(axis)
676               else
677                  tmp_jj = this%boxes(endpoints(j)%box_id)%maxendpoint_id(axis)
678               endif
679     
680               temp = endpoints(i)
681               endpoints(i) = endpoints(j)
682               endpoints(j) = temp
683     
684               if (endpoints(i)%box_id < 0) then
685                  this%boxes(-endpoints(i)%box_id)%minendpoint_id(axis) = tmp_ii
686               else
687                  this%boxes(endpoints(i)%box_id)%maxendpoint_id(axis) = tmp_ii
688               endif
689               if (endpoints(j)%box_id < 0) then
690                  this%boxes(-endpoints(j)%box_id)%minendpoint_id(axis) = tmp_jj
691               else
692                  this%boxes(endpoints(j)%box_id)%maxendpoint_id(axis) = tmp_jj
693               endif
694     
695            elseif (i == j) then
696               marker = i+1
697               return
698            else
699               marker = i
700               return
701            endif
702         end do
703     
704       end subroutine partition
705     
706     end module sweep_and_prune
707