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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module: multi_sweep_and_prune                                       !
4     !                                                                      !
5     !  Purpose: Divides space into a grid of sap_t instances and adds      !
6     !           AABB's to the corresponding sap_t instance(s)              !
7     !                                                                      !
8     !  Reference: http://www.codercorner.com/SAP.pdf                       !
9     !                                                                      !
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
11     
12     #include "version.inc"
13     
14     module multi_sweep_and_prune
15     
16       use sweep_and_prune
17     
18       type multisap_t
19          ! grid particle, e.g. 20x20x20
20          integer :: grid(3)
21     
22          ! bounds of the grid (grid saps on the boundary extend to infinity)
23          real :: minbounds(3), maxbounds(3)
24     
25          ! grid(:)/(maxbounds(:)-minbounds(:))
26          real :: one_over_cell_length(3)
27     
28          ! list saps of length grid(1)*grid(2)*grid(3)
29          type(sap_t), dimension(:), allocatable :: saps
30     
31          ! union of all hashtables for all saps in this multisap
32          type(hashtable_t) :: hashtable
33       end type multisap_t
34     
35       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
36       !                                                                      !
37       !  Type: boxhandle_t                                                   !
38       !                                                                      !
39       !  Purpose: Represents a reference to a particle box in a particular   !
40       !           sap for this multisap.                                     !
41       !                                                                      !
42       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
43     
44       type boxhandle_t
45          integer :: sap_id
46          integer :: box_id
47       end type boxhandle_t
48     
49       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
50       !                                                                      !
51       !  Type: boxhandlelist_t                                               !
52       !                                                                      !
53       !  Purpose: List of boxhandle_t instances that corresponding to a      !
54       !           particular particle with id particle_id.                   !
55       !                                                                      !
56       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
57     
58       ! The maximum number of saps an aabb can belong to. For 2d this could be 4,
59       ! but no harm in just setting it to 8.
60       integer, parameter :: MAX_SAPS = 8
61     
62       type boxhandlelist_t
63          integer :: particle_id
64          type(boxhandle_t) :: list(MAX_SAPS)
65       end type boxhandlelist_t
66     
67      public :: init_multisap, multisap_add, multisap_del, multisap_update, multisap_sort, multisap_quicksort, multisap_sweep, boxhandle_grow
68       private :: multisap_raster
69     
70     contains
71       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
72       !                                                                      !
73       !  Subroutine: init_sap                                                !
74       !                                                                      !
75       !  Purpose: multisap_t constructor                                     !
76       !                                                                      !
77       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
78     
79       subroutine init_multisap(this,x_grid,y_grid,z_grid2,minbounds,maxbounds)
80         use geometry
81         implicit none
82         type(multisap_t), intent(inout) :: this
83         integer, intent(in) :: x_grid,y_grid, z_grid2
84         real, dimension(3), intent(in) :: minbounds, maxbounds
85         integer :: ii,jj,kk, z_grid, sap_id
86     
87         if (NO_K) then
88            z_grid = 1
89         else
90            z_grid = z_grid2
91         endif
92     
93         this%grid(1) = max(x_grid,1)
94         this%grid(2) = max(y_grid,1)
95         this%grid(3) = max(z_grid,1)
96     
97         allocate(this%saps(0:this%grid(1)*this%grid(2)*this%grid(3)-1))
98         do ii=0,this%grid(1)-1
99            do jj=0,this%grid(2)-1
100               do kk=0,this%grid(3)-1
101                  sap_id = (ii*this%grid(2)+jj)*this%grid(3)+kk
102                  call init_sap(this%saps(sap_id),sap_id)
103               enddo
104            enddo
105         enddo
106     
107         this%one_over_cell_length(:) = this%grid(:)/(maxbounds(:)-minbounds(:))
108     
109         this%minbounds(:) = minbounds(:)
110         this%maxbounds(:) = maxbounds(:)
111     
112       end subroutine init_multisap
113     
114       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
115       !                                                                      !
116       !  Subroutine: multisap_raster                                         !
117       !                                                                      !
118       !  Purpose: For a given aabb, return the list sap_ids of ids of the    !
119       !           saps that the aabb intersects according to the grid.       !
120       !           Used by multisap_add and multisap_del.                     !
121       !                                                                      !
122       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
123     
124       subroutine multisap_raster(this,aabb,sap_ids)
125         implicit none
126         type(multisap_t), intent(inout) :: this
127         type(aabb_t), intent(in) :: aabb
128         integer, intent(out) :: sap_ids(MAX_SAPS)
129         integer :: max_grid(3),min_grid(3)
130     
131         sap_ids(:) = -1
132     
133         ! bin to grid
134         min_grid(:) = floor((aabb%minendpoint(:)-this%minbounds(:))*this%one_over_cell_length(:))
135         max_grid(:) = floor((aabb%maxendpoint(:)-this%minbounds(:))*this%one_over_cell_length(:))
136     
137         min_grid(1) = min(max(min_grid(1),0),this%grid(1)-1)
138         min_grid(2) = min(max(min_grid(2),0),this%grid(2)-1)
139         min_grid(3) = min(max(min_grid(3),0),this%grid(3)-1)
140         max_grid(1) = min(max(max_grid(1),0),this%grid(1)-1)
141         max_grid(2) = min(max(max_grid(2),0),this%grid(2)-1)
142         max_grid(3) = min(max(max_grid(3),0),this%grid(3)-1)
143     
144         call add_to_set((min_grid(1)*this%grid(2) + min_grid(2))*this%grid(3) + min_grid(3))
145         call add_to_set((min_grid(1)*this%grid(2) + min_grid(2))*this%grid(3) + max_grid(3))
146         call add_to_set((min_grid(1)*this%grid(2) + max_grid(2))*this%grid(3) + min_grid(3))
147         call add_to_set((min_grid(1)*this%grid(2) + max_grid(2))*this%grid(3) + max_grid(3))
148         call add_to_set((max_grid(1)*this%grid(2) + min_grid(2))*this%grid(3) + min_grid(3))
149         call add_to_set((max_grid(1)*this%grid(2) + min_grid(2))*this%grid(3) + max_grid(3))
150         call add_to_set((max_grid(1)*this%grid(2) + max_grid(2))*this%grid(3) + min_grid(3))
151         call add_to_set((max_grid(1)*this%grid(2) + max_grid(2))*this%grid(3) + max_grid(3))
152     
153       contains
154     
155         !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
156         !                                                                      !
157         !  Subroutine: add_to_set                                              !
158         !                                                                      !
159         !  Purpose: Add sap_id to sap_ids if it's not already in there.        !
160         !                                                                      !
161         !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
162     
163         subroutine add_to_set(sap_id)
164           integer, intent(in) :: sap_id
165           integer :: mm
166     
167           do mm=1, size(sap_ids)
168              if (sap_ids(mm).eq.sap_id) then
169                 ! already in the list
170                 return
171              endif
172              if (sap_ids(mm).eq.-1) then
173                 sap_ids(mm) = sap_id
174                 return
175              endif
176           enddo
177     
178         end subroutine add_to_set
179     
180       end subroutine multisap_raster
181     
182       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
183       !                                                                      !
184       !  Subroutine: multisap_add                                            !
185       !                                                                      !
186       !  Purpose: Add boxes in appropriate saps for aabb representing        !
187       !           particle_id and return list of (sap_id,box_id) handles     !
188       !           in handlelist.                                             !
189       !                                                                      !
190       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
191     
192       subroutine multisap_add(this,aabb,particle_id,handlelist)
193         implicit none
194         type(multisap_t), intent(inout) :: this
195         type(aabb_t), intent(in) :: aabb
196         integer, intent(in) :: particle_id
197         type(boxhandlelist_t), intent(out) :: handlelist
198         integer :: sap_ids(MAX_SAPS)
199         integer :: nn
200     
201         call multisap_raster(this,aabb,sap_ids)
202     
203         handlelist%particle_id = particle_id
204     
205         ! add to each individual SAP
206         do nn=1, size(sap_ids)
207            handlelist%list(nn)%sap_id = sap_ids(nn)
208            if (sap_ids(nn) >= 0) then
209               call add_box(this%saps(sap_ids(nn)),aabb,particle_id,handlelist%list(nn)%box_id)
210            endif
211         enddo
212     
213       end subroutine multisap_add
214     
215       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
216       !                                                                      !
217       !  Subroutine: multisap_del                                            !
218       !                                                                      !
219       !  Purpose: Delete all sap entries corresponding to handlelist.        !
220       !                                                                      !
221       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
222     
223       subroutine multisap_del(this,handlelist)
224         implicit none
225         type(multisap_t), intent(inout) :: this
226         type(boxhandlelist_t), intent(in) :: handlelist
227         integer :: nn
228     
229         ! remove from each SAP listed for this id
230         do nn=1, size(handlelist%list)
231            if (handlelist%list(nn)%sap_id >= 0) then
232               call del_box(this%saps(handlelist%list(nn)%sap_id),handlelist%list(nn)%box_id)
233            endif
234         enddo
235     
236       end subroutine multisap_del
237     
238       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
239       !                                                                      !
240       !  Subroutine: multisap_update                                         !
241       !                                                                      !
242       !  Purpose: Update the sap entries corresponding to handlelist with    !
243       !           new location aabb.  Returns an updated handlelist with     !
244       !           possibly different sap entries.                            !
245       !                                                                      !
246       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
247     
248       subroutine multisap_update(this,aabb,handlelist)
249         implicit none
250         type(multisap_t), intent(inout) :: this
251         type(aabb_t), intent(in) :: aabb
252         type(boxhandlelist_t), intent(inout) :: handlelist
253         integer, DIMENSION(MAX_SAPS) :: new_sap_ids
254         logical :: found
255         integer :: mm,nn,first_blank, box_id
256     
257         call multisap_raster(this,aabb,new_sap_ids)
258     
259         ! update for each SAP listed for this id
260         do nn=1, size(new_sap_ids)
261            if (new_sap_ids(nn) < 0) exit
262     
263            found = .false.
264            first_blank = -1
265            do mm=1, size(handlelist%list)
266               if (handlelist%list(mm)%sap_id < 0) first_blank = mm
267     
268               if (handlelist%list(mm)%sap_id .eq. new_sap_ids(nn)) then
269                  ! update existing SAPs
270                  box_id = handlelist%list(mm)%box_id
271                  call update_box(this%saps(new_sap_ids(nn)),handlelist%list(mm)%box_id,aabb)
272                  found = .true.
273                  exit
274               endif
275            enddo
276            if (.not.found) then
277               ! add SAPs yet not listed for this id
278     
279               if (first_blank .eq. -1) then
280                  print *,"AABB is in more than 8 saps? Something is wrong:  ",first_blank,handlelist%list
281                  ERROR_STOP __LINE__
282               endif
283               call add_box(this%saps(new_sap_ids(nn)),aabb,handlelist%particle_id,handlelist%list(first_blank)%box_id)
284               handlelist%list(first_blank)%sap_id = new_sap_ids(nn)
285            endif
286         enddo
287     
288         do mm=1, size(handlelist%list)
289            if (handlelist%list(mm)%sap_id < 0) cycle
290     
291            found = .false.
292            do nn=1, size(new_sap_ids)
293               if (new_sap_ids(nn) < 0) exit
294     
295               if (handlelist%list(mm)%sap_id .eq. new_sap_ids(nn)) then
296                  found = .true.
297                  exit
298               endif
299            enddo
300     
301            if (.not.found) then
302               ! remove SAP not found in handle%sap_id
303               call del_box(this%saps(handlelist%list(mm)%sap_id),handlelist%list(mm)%box_id)
304               handlelist%list(mm)%sap_id = -1
305            endif
306     
307         enddo
308     
309       end subroutine multisap_update
310     
311       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
312       !                                                                      !
313       !  Subroutine: multisap_sort                                           !
314       !                                                                      !
315       !  Purpose: Call sort on each of the saps in this multisap,            !
316       !           then set the multisap's hashtable to the union of the      !
317       !           hashtables of all the saps.                                !
318       !                                                                      !
319       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
320     
321       subroutine multisap_sort(this)
322         ! use discretelement
323         ! use functions, only: is_nonexistent
324         implicit none
325         type(multisap_t), intent(inout) :: this
326         integer :: ii, jj, kk, sap_id
327         integer :: pair(2)
328         integer :: sighs, cc,ll,i, cc_start, cc_end
329     
330         call init_pairs(this%hashtable)
331     
332     !$omp parallel default(none) private(ii,jj,kk,sap_id,pair) shared(this)
333     !$omp do collapse(3)
334         do ii=0,this%grid(1)-1
335            do jj=0,this%grid(2)-1
336               do kk=0,this%grid(3)-1
337                  sap_id = ii*this%grid(2)*this%grid(3)+jj*this%grid(3)+kk
338                  call sort(this%saps(sap_id))
339               enddo
340            enddo
341         enddo
342     !$omp end parallel
343     
344         print *," NUMBER OF DESGRIDS IS: ",this%grid(1)*this%grid(2)*this%grid(3)
345     
346      sighs = 0
347     
348      open (unit=235,file="pairs.txt",action="write",status="replace")
349     
350        do ii=0,this%grid(1)-1
351            do jj=0,this%grid(2)-1
352               do kk=0,this%grid(3)-1
353                  sap_id = ii*this%grid(2)*this%grid(3)+jj*this%grid(3)+kk
354                  call reset_pairs(this%saps(sap_id)%hashtable)
355                  sighs = sighs + this%saps(sap_id)%hashtable%table_size
356                  do
357                     call get_pair(this%saps(sap_id)%hashtable,pair)
358                     if (pair(1).eq.0 .and. pair(2).eq.0) exit
359                     write (235,*) pair(1),pair(2)
360                     call add_pair(this%hashtable,pair(1),pair(2))
361     
362                  enddo
363               enddo
364            enddo
365         enddo
366     
367         close (unit=235)
368     
369         ! call reset_pairs(this%hashtable)
370     
371         ! do
372         !    call get_pair(this%hashtable,pair)
373         !    if (pair(1).eq.0 .and. pair(2).eq.0) exit
374         !    write (235,*) pair(1),pair(2)
375         ! enddo
376     
377         ! print *,"multisize ",this%hashtable%table_size
378         ! print *,"sighs ",sighs
379         ! stop __LINE__
380     
381       end subroutine multisap_sort
382     
383       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
384       !                                                                      !
385       !  Subroutine: multisap_add                                            !
386       !                                                                      !
387       !  Purpose: Call quicksort on each of the saps in this multisap.       !
388       !                                                                      !
389       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
390     
391       subroutine multisap_quicksort(this)
392         implicit none
393         type(multisap_t), intent(inout) :: this
394         integer :: ii, jj, kk
395         integer :: sap_id, boxcount
396     
397         boxcount = 0
398     
399         do ii=0,this%grid(1)-1
400            do jj=0,this%grid(2)-1
401               do kk=0,this%grid(3)-1
402                  sap_id = ii*this%grid(2)*this%grid(3)+jj*this%grid(3)+kk
403                  call quicksort(this%saps(sap_id))
404                  boxcount = boxcount + this%saps(sap_id)%boxes_len
405               enddo
406            enddo
407         enddo
408     
409       end subroutine multisap_quicksort
410     
411       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
412       !                                                                      !
413       !  Subroutine: multisap_sweep                                          !
414       !                                                                      !
415       !  Purpose: Call sweep on each of the saps in this multisap.           !
416       !                                                                      !
417       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
418     
419       subroutine multisap_sweep(this)
420         implicit none
421         type(multisap_t), intent(inout) :: this
422         integer :: ii, jj, kk
423         integer :: sap_id
424     
425         do ii=0,this%grid(1)-1
426            do jj=0,this%grid(2)-1
427               do kk=0,this%grid(3)-1
428                  sap_id = (ii*this%grid(2)+jj)*this%grid(3)+kk
429                  call sweep(this%saps(sap_id))
430               enddo
431            enddo
432         enddo
433     
434       end subroutine multisap_sweep
435     
436       subroutine multisap_check(this)
437         implicit none
438         type(multisap_t), intent(inout) :: this
439         integer :: ii, jj, kk
440         integer :: sap_id
441     
442         do ii=0,this%grid(1)-1
443            do jj=0,this%grid(2)-1
444               do kk=0,this%grid(3)-1
445                  sap_id = (ii*this%grid(2)+jj)*this%grid(3)+kk
446     
447                  call check(this%saps(sap_id))
448               enddo
449            enddo
450         enddo
451     
452       end subroutine multisap_check
453     
454       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
455       !                                                                      !
456       !  Subroutine: boxhandle_GROW                                          !
457       !                                                                      !
458       !  Purpose: resize array of box_t                                      !
459       !                                                                      !
460       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
461     
462       subroutine boxhandle_GROW(boxhandles,new_size)
463         IMPLICIT NONE
464     
465         INTEGER, INTENT(IN) :: new_size
466         type(boxhandlelist_t), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: boxhandles
467         type(boxhandlelist_t), DIMENSION(:), ALLOCATABLE :: boxhandle_tmp
468         INTEGER lSIZE
469     
470         lSIZE = size(boxhandles,1)
471         allocate(boxhandle_tmp(new_size))
472         boxhandle_tmp(1:lSIZE) = boxhandles(1:lSIZE)
473         call move_alloc(boxhandle_tmp,boxhandles)
474     
475       end subroutine boxhandle_GROW
476     
477     end module multi_sweep_and_prune
478