File: /nfs/home/0/users/jenkins/mfix.git/model/des/layout_mi_dem.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: SET_BC_DEM                                              !
4     !  Author: J.Musser                                   Date: 13-Jul-09  !
5     !                                                                      !
6     !  Purpose: Check the data provided for the des mass inflow boundary   !
7     !  condition and flag errors if the data is improper.  This module is  !
8     !  also used to convert the proveded information into the format       !
9     !  necessary for the dependent subrountines to function properly.      !
10     !                                                                      !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
12           SUBROUTINE LAYOUT_MI_DEM(BCV, BCV_I, MAX_DIA)
13     
14           use bc, only: BC_PLANE
15           USE run, only: RUN_TYPE
16     
17           use error_manager
18     
19           IMPLICIT NONE
20     
21           INTEGER, INTENT(IN) :: BCV
22           INTEGER, INTENT(IN) :: BCV_I      ! BC loop counter
23           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
24     ! Local debug flags.
25           LOGICAL, parameter :: setDBG = .FALSE.
26           LOGICAL, parameter :: showMAP = .FALSE.
27     
28           CALL INIT_ERR_MSG("LAYOUT_MI_DEM")
29     
30     ! This subroutine determines the pattern that the particles will need to
31     ! enter the system, if any. This routine only needs to be called if a
32     ! run is new.  If a run is a RESTART_1, all of the setup information
33     ! provided by this subroutine is will be obtained from the *_DES.RES file.
34     ! This is done due to this routine's strong dependence on the
35     ! RANDOM_NUMBER() subroutine.
36           IF(RUN_TYPE == 'RESTART_1') THEN
37              CALL SET_DEM_MI_OWNER(BCV, BCV_I)
38           ELSE
39              SELECT CASE (BC_PLANE(BCV))
40              CASE('N','S')
41                 CALL LAYOUT_DEM_MI_NS(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
42              CASE('E','W')
43                 CALL LAYOUT_DEM_MI_EW(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
44              CASE('T','B')
45                 CALL LAYOUT_DEM_MI_TB(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
46              END SELECT
47           ENDIF
48     
49           CALL FINL_ERR_MSG
50     
51           RETURN
52           END SUBROUTINE LAYOUT_MI_DEM
53     
54     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
55     !                                                                      !
56     !  Subroutine: LAYOUT_DEM_MI_NS                                        !
57     !                                                                      !
58     !  Purpose:  This routine determines the layout of the mass inlet as   !
59     !  either ordered or random based upon the inlet conditions.  This     !
60     !  routine also verifies that the specified inlet conditions for mass  !
61     !  or volumetric flow rates along with inlet size (length or area) and !
62     !  particle inlet velocity will work.  If not an error is flagged and  !
63     !  the program is exited.                                              !
64     !                                                                      !
65     !                                                                      !
66     !  Author: J.Musser                                   Date: 14-Aug-09  !
67     !                                                                      !
68     !  Comments:                                                           !
69     !                                                                      !
70     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
71           SUBROUTINE LAYOUT_DEM_MI_NS(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
72     
73     ! Global variables:
74     !---------------------------------------------------------------------//
75           use bc, only: BC_PLANE, BC_Y_s, BC_AREA
76           use bc, only: BC_X_w, BC_X_e
77           use bc, only: BC_Z_b, BC_Z_t
78     
79           use des_bc, only: DEM_MI
80     
81           use stl, only: N_FACETS_DES
82           use stl, only: VERTEX, NORM_FACE
83           use cutcell, only: USE_STL
84     
85           use compar, only: myPE
86           use geometry, only: IMAX, JMAX, KMAX
87           use geometry, only: DX, DY, DZ
88           use geometry, only: XMIN, DO_K
89     
90           use funits, only: DMP_LOG
91     
92     ! Global parameters:
93     !---------------------------------------------------------------------//
94           use param1, only: ZERO, HALF, ONE
95     
96     ! Module procedures
97     !---------------------------------------------------------------------//
98           use des_bc, only: EXCLUDE_DEM_MI_CELL
99           use mpi_utility, only: GLOBAL_ALL_SUM
100           use des_stl_functions, only: TestTriangleAABB
101           use functions, only: IS_ON_myPE_OWNS
102     
103           use error_manager
104     
105           IMPLICIT NONE
106     
107     ! Dummy arguments
108     !---------------------------------------------------------------------//
109     ! passed arguments giving index information of the boundary
110           INTEGER, INTENT(IN) :: BCV
111           INTEGER, INTENT(IN) :: BCV_I
112     ! Max diameter of incoming particles at bc
113           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
114     ! Debug flas.
115           LOGICAL, INTENT(IN) :: setDBG, showMAP
116     
117     ! Local variables
118     !---------------------------------------------------------------------//
119     ! Loop counters.
120           INTEGER LL, LC
121     ! Indices for mapping to fluid grid.
122           INTEGER I, J, K
123     ! Local MESH indices
124           INTEGER H, W
125     ! Temp variable for double precision
126           DOUBLE PRECISION :: TMP_DP
127     ! Temporary variable for integers
128           INTEGER :: TMP_INT
129     ! Window indicies.
130           INTEGER, allocatable :: MESH_H(:)
131           INTEGER, allocatable :: MESH_W(:)
132     ! Window positions.
133           DOUBLE PRECISION, allocatable :: MESH_P(:)
134           DOUBLE PRECISION, allocatable :: MESH_Q(:)
135     ! Map of array index to MI cell
136           INTEGER, allocatable :: RAND_MAP(:)
137     ! Map of MI cells to owner processes
138           INTEGER, allocatable :: FULL_MAP(:,:)
139     ! max number of partitions along length of inlet
140           INTEGER :: WMAX, HMAX
141     ! the length of each side of the inlet boundary
142           DOUBLE PRECISION :: PLEN, QLEN
143     ! Number of occupied mesh cells
144           INTEGER :: OCCUPANTS
145     ! Offset and window size.
146           DOUBLE PRECISION :: SHIFT, WINDOW
147     ! The origin and dimenion of MI cells. (STL intersection tests)
148           DOUBLE PRECISION :: ORIGIN(3), EXTENTS(3)
149     ! Separating axis test dummy variable
150           INTEGER :: SEP_AXIS
151     ! Indicates that a separating axis exists
152           LOGICAL :: SA_EXIST
153     ! Debug flag.
154           LOGICAL :: dFlag
155     !......................................................................!
156     
157     ! Initialize the error manager.
158           CALL INIT_ERR_MSG('LAYOUT_DEM_MI_NS')
159     
160           dFlag = (DMP_LOG .AND. setDBG)
161           if(dFlag) write(*,"(2/,'Building NS DEM_MI: ',I3)") BCV_I
162     
163     ! Store the index that maps back to the user input.
164     
165           OCCUPANTS = 0
166     
167     ! Calculate number of partitions in first direction.
168           PLEN = BC_X_e(BCV) - BC_X_w(BCV)
169           WMAX = FLOOR(real(PLEN/MAX_DIA))
170           allocate( MESH_W(WMAX) )
171           allocate( MESH_P(0:WMAX) )
172     
173     ! Calculate number of partitions in the second direction.
174           QLEN = merge(BC_Z_t(BCV) - BC_Z_b(BCV), MAX_DIA, DO_K)
175           HMAX = FLOOR(real(QLEN/MAX_DIA))
176           allocate( MESH_H(HMAX) )
177           allocate( MESH_Q(0:HMAX) )
178     
179     ! Allocate the full map.
180           allocate( FULL_MAP(WMAX, HMAX))
181     
182     ! Set the value of the boundary condtion offset value used in the
183     ! placement of new particles.
184           CALL CALC_CELL_INTERSECT(ZERO, BC_Y_s(BCV), DY, JMAX, J)
185           SHIFT = merge(-ONE, ONE, BC_PLANE(BCV) == 'N')
186           DEM_MI(BCV_I)%OFFSET = BC_Y_s(BCV) + MAX_DIA*SHIFT
187           DEM_MI(BCV_I)%L = J + int(SHIFT)
188           if(dFlag) write(*,"(2x,'Offset: ',3x,I4,3x,g12.5)") &
189              DEM_MI(BCV_I)%L, DEM_MI(BCV_I)%OFFSET
190     
191     
192     ! Dimension of grid cell for seeding particles; this may be larger than
193     ! than the particle diameter but not smaller:
194           DEM_MI(BCV_I)%WINDOW = MIN(PLEN/WMAX, QLEN/HMAX)
195           WINDOW = DEM_MI(BCV_I)%WINDOW
196           if(dFlag) write(*,"(2x,'Windows size: ',g12.5)") WINDOW
197     
198     ! Setup the first direction.
199           SHIFT = HALF*(PLEN - WMAX*WINDOW)
200           MESH_P(0) = BC_X_w(BCV) + SHIFT
201           if(dFlag) write(*,8005) 'P', SHIFT, 'P', MESH_P(0)
202           DO LC=1,WMAX
203              MESH_P(LC) = MESH_P(0) + dble(LC-1)*WINDOW
204              SHIFT = MESH_P(LC) + HALF*WINDOW
205              CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, MESH_W(LC))
206              IF(dFlag)WRITE(*,8006) LC, 'W', MESH_W(LC), 'P', MESH_P(LC)
207           ENDDO
208     
209     ! Setup the second direction.
210           IF(DO_K) THEN
211              SHIFT = HALF*(QLEN - HMAX*WINDOW)
212              MESH_Q(0) = BC_Z_b(BCV) + SHIFT
213              if(dFlag) write(*,8005) 'Q',SHIFT, 'Q',MESH_Q(0)
214              DO LC=1,HMAX
215                 MESH_Q(LC) = MESH_Q(0) + dble(LC-1)*WINDOW
216                 SHIFT = MESH_Q(LC) + HALF*WINDOW
217                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, MESH_H(LC))
218                 IF(dFlag)WRITE(*,8006) LC, 'H', MESH_H(LC), 'Q', MESH_Q(LC)
219              ENDDO
220           ELSE
221              MESH_H(1) = 1
222              MESH_Q(1) = 0.0d0
223           ENDIF
224     
225     ! Get the Jth index of the fluid cell
226           CALL CALC_CELL_INTERSECT(ZERO, BC_Y_s(BCV), DY, JMAX, J)
227     
228     
229     ! If the computationsl cell adjacent to the DEM_MI mesh cell is a
230     ! fluid cell and has not been cut, store the ID of the cell owner.
231           DO H=1,HMAX
232           DO W=1,WMAX
233     
234              I = MESH_W(W)
235              K = MESH_H(H)
236              FULL_MAP(W,H) = 0
237     
238              IF(.NOT.IS_ON_myPE_owns(I,J,K)) CYCLE
239     
240              IF(DO_K) THEN
241     
242                 CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
243                 CALL CALC_CELL_INTERSECT(ZERO, MESH_Q(H), DZ, KMAX, K)
244                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
245     
246                 SHIFT = MESH_P(W)+WINDOW
247                 CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, I)
248                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
249     
250                 SHIFT = MESH_Q(H)+WINDOW
251                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, K)
252                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
253     
254                 CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
255                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
256     
257              ELSE
258     
259                 K = MESH_H(1)
260                 CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
261                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
262     
263                 SHIFT = MESH_P(W)+WINDOW
264                 CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, I)
265                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
266              ENDIF
267     
268              FULL_MAP(W,H) = myPE+1
269           ENDDO
270           ENDDO
271     
272     
273     ! For complex boundaries defined by STLs, exclude any mass inflow cells
274     ! that intersect the boundary and all cells opposite the normal. The MI
275     ! cell sizes are increased by 10% to provide a small buffer.
276           IF(USE_STL) THEN
277     
278              EXTENTS(2) = 6.0d0*MAX_DIA
279              EXTENTS(1) = WINDOW * 1.10d0
280              EXTENTS(3) = WINDOW * 1.10d0
281     
282              DO H=1,HMAX
283              DO W=1,WMAX
284     
285                 ORIGIN(2) = BC_Y_s(BCV) - 3.0d0*MAX_DIA
286                 ORIGIN(1) = MESH_P(W) - WINDOW * 0.05d0
287                 ORIGIN(3) = MESH_Q(H) - WINDOW * 0.05d0
288     
289                 FACET_LP: DO LC=1, N_FACETS_DES
290     
291                    IF(BC_Y_s(BCV) > maxval(VERTEX(:,2,LC))) CYCLE FACET_LP
292                    IF(BC_Y_s(BCV) < minval(VERTEX(:,2,LC))) CYCLE FACET_LP
293     
294                    IF(BC_X_w(BCV) > maxval(VERTEX(:,1,LC))) CYCLE FACET_LP
295                    IF(BC_X_e(BCV) < minval(VERTEX(:,1,LC))) CYCLE FACET_LP
296     
297                    IF(BC_Z_b(BCV) > maxval(VERTEX(:,3,LC))) CYCLE FACET_LP
298                    IF(BC_Z_t(BCV) < minval(VERTEX(:,3,LC))) CYCLE FACET_LP
299     
300                    CALL TESTTRIANGLEAABB(VERTEX(:,:,LC), NORM_FACE(:,LC),  &
301                       ORIGIN(:), EXTENTS(:), SA_EXIST, SEP_AXIS, I, J, K)
302     
303                    IF(.NOT.SA_EXIST) THEN
304                       IF(NORM_FACE(1,LC) >= 0) THEN
305                          FULL_MAP(1:W,H) = 0
306                       ELSE
307                          FULL_MAP(W:WMAX,H) = 0
308                       ENDIF
309                       IF(NORM_FACE(3,LC) >= 0) THEN
310                          FULL_MAP(W,1:H) = 0
311                       ELSE
312                          FULL_MAP(W,H:HMAX) = 0
313                       ENDIF
314                    ENDIF
315                 ENDDO FACET_LP
316              ENDDO
317              ENDDO
318           ENDIF
319     
320     ! Add up the total number of available positions in the seeding map.
321           DO H=1,HMAX
322           DO W=1,WMAX
323              IF(FULL_MAP(W,H) /= 0) OCCUPANTS = OCCUPANTS + 1
324           ENDDO
325           ENDDO
326     
327     ! Sync the full map across all ranks.
328           CALL GLOBAL_ALL_SUM(OCCUPANTS)
329           CALL GLOBAL_ALL_SUM(FULL_MAP)
330     
331     ! Throw an error and exit if there are no occupants.
332           IF(OCCUPANTS == 0) THEN
333              WRITE(ERR_MSG, 1100) BCV_I
334              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
335           ENDIF
336     
337      1100 FORMAT('Error 1100: No un-cut fluid cells adjacent to DEM_MI ',  &
338              'staging area.',/'Unable to setup the discrete solids mass ', &
339              'inlet for BC:',I3)
340     
341     ! Store the number of occupants.
342           DEM_MI(BCV_I)%OCCUPANTS = OCCUPANTS
343     
344     ! Display the fill map if debugging
345           IF(dFlag .OR. showMAP) THEN
346              WRITE(*,"(2/,2x,'Displaying Fill Map:')")
347              DO H=HMAX,1,-1
348                 WRITE(*,"(2x,'H =',I3)",advance='no')H
349                 DO W=1,WMAX
350                    IF(FULL_MAP(W,H) == 0) then
351                       WRITE(*,"(' *')",advance='no')
352                    ELSE
353                       WRITE(*,"(' .')",advance='no')
354                    ENDIF
355                 ENDDO
356                 WRITE(*,*)' '
357              ENDDO
358           ENDIF
359     
360     ! Construct an array of integers randomly ordered.
361           if(dFLAG) write(*,"(2/,2x,'Building RAND_MAP:')")
362           allocate( RAND_MAP(OCCUPANTS) )
363           RAND_MAP = 0
364     
365     ! Only Rank 0 will randomize the layout and broadcast it to the other
366     ! ranks. This will ensure that all ranks see the same layout.
367           IF(myPE == 0) THEN
368              LL = 1
369              DO WHILE (RAND_MAP(OCCUPANTS) .EQ. 0)
370                 CALL RANDOM_NUMBER(TMP_DP)
371                 TMP_INT = CEILING(real(TMP_DP*dble(OCCUPANTS)))
372                 DO LC = 1, LL
373                   IF(TMP_INT .EQ. RAND_MAP(LC) )EXIT
374                   IF(LC .EQ. LL)THEN
375                      if(dFlag) WRITE(*,"(4x,'LC:',I9,' : ',I9)") LC, TMP_INT
376                      RAND_MAP(LC) = TMP_INT
377                      LL = LL + 1
378                   ENDIF
379                 ENDDO
380              ENDDO
381           ENDIF
382     
383           CALL GLOBAL_ALL_SUM(RAND_MAP)
384     
385     ! Initialize the vacancy pointer.
386           DEM_MI(BCV_I)%VACANCY = 1
387     
388     ! Allocate the mass inlet storage arrays.
389           allocate( DEM_MI(BCV_I)%W(OCCUPANTS) )
390           allocate( DEM_MI(BCV_I)%P(OCCUPANTS) )
391           allocate( DEM_MI(BCV_I)%H(OCCUPANTS) )
392           allocate( DEM_MI(BCV_I)%Q(OCCUPANTS) )
393           allocate( DEM_MI(BCV_I)%OWNER(OCCUPANTS) )
394     
395           if(dFlag) write(*,8010)
396     ! Store the MI layout in the randomized order.
397           LC = 0
398           DO H=1,HMAX
399           DO W=1,WMAX
400              IF(FULL_MAP(W,H) == 0) CYCLE
401              LC = LC + 1
402              LL = RAND_MAP(LC)
403              DEM_MI(BCV_I)%OWNER(LL) = FULL_MAP(W,H) - 1
404     
405              DEM_MI(BCV_I)%W(LL) = MESH_W(W)
406              DEM_MI(BCV_I)%H(LL) = MESH_H(H)
407     
408              DEM_MI(BCV_I)%P(LL) = MESH_P(W)
409              DEM_MI(BCV_I)%Q(LL) = MESH_Q(H)
410     
411              if(dFlag) write(*,8011) DEM_MI(BCV_I)%OWNER(LL), &
412                 DEM_MI(BCV_I)%W(LL), DEM_MI(BCV_I)%H(LL), DEM_MI(BCV_I)%L, &
413                 DEM_MI(BCV_I)%P(LL), DEM_MI(BCV_I)%Q(LL), DEM_MI(BCV_I)%OFFSET
414     
415           ENDDO
416           ENDDO
417     
418     
419      8010 FORMAT(2/,2x,'Storing DEM_MI data:',/4X,'OWNER',5X,'W',5X,'H',   &
420              5X,'L',7X,'P',12X,'Q',12X,'R')
421      8011 FORMAT(4x,I5,3(2X,I4),3(2x,g12.5))
422     
423     
424           if(dFlag) write(*,"(2/,2x,'Inlet area sizes:')")
425           if(dFlag) write(*,9000) 'mfix.dat: ', PLEN * QLEN
426           if(dFlag) write(*,9000) 'BC_AREA:  ', BC_AREA(BCV)
427           if(dFlag) write(*,9000) 'DEM_MI:   ', OCCUPANTS * (WINDOW**2)
428      9000 FORMAT(2x,A,g12.5)
429     
430     ! House keeping.
431           IF( allocated(MESH_H)) deallocate(MESH_H)
432           IF( allocated(MESH_W)) deallocate(MESH_W)
433           IF( allocated(MESH_P)) deallocate(MESH_P)
434           IF( allocated(MESH_Q)) deallocate(MESH_Q)
435     
436           IF( allocated(RAND_MAP)) deallocate(RAND_MAP)
437           IF( allocated(FULL_MAP)) deallocate(FULL_MAP)
438     
439           CALL FINL_ERR_MSG
440     
441           RETURN
442     
443      8005 FORMAT(2/,2x,'Building MESH_',A1,':',/4x,'Shift:',f8.4,/4x,      &
444              'MESH_',A1,'(0) = ',f8.4,/)
445     
446      8006 FORMAT(4x,'LC = ',I4,3x,A1,' =',I3,3x,A1,' =',f8.4)
447     
448           END SUBROUTINE LAYOUT_DEM_MI_NS
449     
450     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
451     !                                                                      !
452     !  Subroutine: LAYOUT_DEM_MI_EW                                        !
453     !                                                                      !
454     !  Purpose:  This routine determines the layout of the mass inlet as   !
455     !  either ordered or random based upon the inlet conditions.  This     !
456     !  routine also verifies that the specified inlet conditions for mass  !
457     !  or volumetric flow rates along with inlet size (length or area) and !
458     !  particle inlet velocity will work.  If not an error is flagged and  !
459     !  the program is exited.                                              !
460     !                                                                      !
461     !                                                                      !
462     !  Author: J.Musser                                   Date: 14-Aug-09  !
463     !                                                                      !
464     !  Comments:                                                           !
465     !                                                                      !
466     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
467           SUBROUTINE LAYOUT_DEM_MI_EW(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
468     
469           use bc, only: BC_PLANE, BC_X_w, BC_AREA
470           use bc, only: BC_Y_s, BC_Y_n
471           use bc, only: BC_Z_b, BC_Z_t
472     
473           use des_bc, only: DEM_MI
474     
475           use stl, only: N_FACETS_DES
476           use stl, only: VERTEX, NORM_FACE
477           use cutcell, only: USE_STL
478     
479           use compar, only: myPE
480           use geometry, only: IMAX, JMAX, KMAX
481           use geometry, only: DX, DY, DZ
482           use geometry, only: XMIN, DO_K
483     
484           use funits, only: DMP_LOG
485     
486     ! Global parameters:
487     !---------------------------------------------------------------------//
488           use param1, only: ZERO, HALF, ONE
489     
490     ! Module procedures
491     !---------------------------------------------------------------------//
492           use des_bc, only: EXCLUDE_DEM_MI_CELL
493           use des_stl_functions, only: TestTriangleAABB
494           use mpi_utility, only: GLOBAL_ALL_SUM
495           use functions, only: IS_ON_myPE_OWNS
496     
497           use error_manager
498     
499           IMPLICIT NONE
500     
501     ! Dummy arguments
502     !---------------------------------------------------------------------//
503     ! passed arguments giving index information of the boundary
504           INTEGER, INTENT(IN) :: BCV
505           INTEGER, INTENT(IN) :: BCV_I
506     ! Max diameter of incoming particles at bc
507           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
508     ! Debug flas.
509           LOGICAL, INTENT(IN) :: setDBG, showMAP
510     
511     ! Local variables
512     !---------------------------------------------------------------------//
513     ! Loop counters.
514           INTEGER LL, LC
515     ! Indices for mapping to fluid grid.
516           INTEGER IJK, I, J, K
517     ! Local MESH indices
518           INTEGER H, W
519     ! Temp variable for double precision
520           DOUBLE PRECISION :: TMP_DP
521     ! Temporary variable for integers
522           INTEGER :: TMP_INT
523     ! Window indicies.
524           INTEGER, allocatable :: MESH_H(:)
525           INTEGER, allocatable :: MESH_W(:)
526     ! Window positions.
527           DOUBLE PRECISION, allocatable :: MESH_P(:)
528           DOUBLE PRECISION, allocatable :: MESH_Q(:)
529     ! Map of array index to MI cell
530           INTEGER, allocatable :: RAND_MAP(:)
531     ! Map of MI cells to owner processes
532           INTEGER, allocatable :: FULL_MAP(:,:)
533     ! max number of partitions along length of inlet
534           INTEGER :: WMAX, HMAX
535     ! the length of each side of the inlet boundary
536           DOUBLE PRECISION :: PLEN, QLEN
537     ! Number of occupied mesh cells
538           INTEGER :: OCCUPANTS
539     ! Offset and window size.
540           DOUBLE PRECISION :: SHIFT, WINDOW
541     ! The origin and dimenion of MI cells. (STL intersection tests)
542           DOUBLE PRECISION :: ORIGIN(3), EXTENTS(3)
543     ! Separating axis test dummy variable
544           INTEGER :: SEP_AXIS
545     ! Indicates that a separating axis exists
546           LOGICAL :: SA_EXIST
547     ! Local debug flag.
548           LOGICAL :: dFlag
549     !......................................................................!
550     
551     
552     ! Initialize the error manager.
553           CALL INIT_ERR_MSG('LAYOUT_DEM_MI_EW')
554     
555           dFlag = (DMP_LOG .AND. setDBG)
556           if(dFlag) write(*,"(2/,'Building EW DEM_MI: ',I3)") BCV_I
557     
558           OCCUPANTS = 0
559     
560     ! Calculate number of partitions in first direction.
561           PLEN = BC_Y_n(BCV) - BC_Y_s(BCV)
562           WMAX = FLOOR(real(PLEN/MAX_DIA))
563           allocate( MESH_W(WMAX) )
564           allocate( MESH_P(0:WMAX) )
565     
566     ! Calculate number of partitions in the second direction.
567           QLEN = merge(BC_Z_t(BCV) - BC_Z_b(BCV), MAX_DIA, DO_K)
568           HMAX = FLOOR(real(QLEN/MAX_DIA))
569           allocate( MESH_H(HMAX) )
570           allocate( MESH_Q(0:HMAX) )
571     
572     ! Allocate the full map.
573           allocate( FULL_MAP(WMAX, HMAX))
574     
575     ! Set the value of the boundary condtion offset value used in the
576     ! placement of new particles.
577           CALL CALC_CELL_INTERSECT(XMIN, BC_X_w(BCV), DX, IMAX, I)
578           SHIFT = merge(-ONE, ONE, BC_PLANE(BCV) == 'E')
579           DEM_MI(BCV_I)%OFFSET = BC_X_w(BCV) + MAX_DIA*SHIFT
580           DEM_MI(BCV_I)%L = I + int(SHIFT)
581           if(dFlag) write(*,"(2x,'Offset: ',3x,I4,3x,g12.5)") &
582              DEM_MI(BCV_I)%L, DEM_MI(BCV_I)%OFFSET
583     
584     
585     ! Dimension of grid cell for seeding particles; this may be larger than
586     ! than the particle diameter but not smaller:
587           DEM_MI(BCV_I)%WINDOW = MIN(PLEN/WMAX, QLEN/HMAX)
588           WINDOW = DEM_MI(BCV_I)%WINDOW
589           if(dFlag) write(*,"(2x,'Windows size: ',g12.5)") WINDOW
590     
591     ! Setup the first direction.
592           SHIFT = HALF*(PLEN - WMAX*WINDOW)
593           MESH_P(0) = BC_Y_s(BCV) + SHIFT
594           if(dFlag) write(*,8005) 'P', SHIFT, 'P', MESH_P(0)
595           DO LC=1,WMAX
596              MESH_P(LC) = MESH_P(0) + dble(LC-1)*WINDOW
597              SHIFT = MESH_P(LC) + HALF*WINDOW
598              CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, MESH_W(LC))
599              IF(dFlag)WRITE(*,8006) LC, 'W', MESH_W(LC), 'P', MESH_P(LC)
600           ENDDO
601     
602     ! Setup the second direction.
603           IF(DO_K) THEN
604              SHIFT = HALF*(QLEN - HMAX*WINDOW)
605              MESH_Q(0) = BC_Z_b(BCV) + SHIFT
606              if(dFlag) write(*,8005) 'Q',SHIFT, 'Q',MESH_Q(0)
607              DO LC=1,HMAX
608                 MESH_Q(LC) = MESH_Q(0) + dble(LC-1)*WINDOW
609                 SHIFT = MESH_Q(LC) + HALF*WINDOW
610                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, MESH_H(LC))
611                 IF(dFlag)WRITE(*,8006) LC, 'H', MESH_H(LC), 'Q', MESH_Q(LC)
612              ENDDO
613           ELSE
614              MESH_H(1) = 1
615              MESH_Q(1) = 0.0d0
616           ENDIF
617     
618     ! Get the Jth index of the fluid cell
619           CALL CALC_CELL_INTERSECT(XMIN, BC_X_w(BCV), DX, IMAX, I)
620     
621     
622     ! If the computationsl cell adjacent to the DEM_MI mesh cell is a
623     ! fluid cell, store the rank of the cell owner.
624           DO H=1,HMAX
625           DO W=1,WMAX
626     
627              J = MESH_W(W)
628              K = MESH_H(H)
629              FULL_MAP(W,H) = 0
630     
631              IF(.NOT.IS_ON_myPE_owns(I,J,K)) CYCLE
632     
633              IF(DO_K) THEN
634     
635                 CALL CALC_CELL_INTERSECT(ZERO, MESH_P(W), DY, JMAX, J)
636                 CALL CALC_CELL_INTERSECT(ZERO, MESH_Q(H), DZ, KMAX, K)
637                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
638     
639                 SHIFT = MESH_P(W)+WINDOW
640                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, J)
641                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
642     
643                 SHIFT = MESH_Q(H)+WINDOW
644                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, K)
645                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
646     
647                 CALL CALC_CELL_INTERSECT(ZERO, MESH_P(W), DY, JMAX, J)
648                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
649     
650              ELSE
651     
652                 K = MESH_H(1)
653                 CALL CALC_CELL_INTERSECT(ZERO, MESH_P(W), DY, JMAX, J)
654                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
655     
656                 SHIFT = MESH_P(W)+WINDOW
657                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, J)
658                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
659              ENDIF
660     
661              FULL_MAP(W,H) = myPE+1
662           ENDDO
663           ENDDO
664     
665     
666     ! For complex boundaries defined by STLs, exclude any mass inflow cells
667     ! that intersect the boundary and all cells opposite the normal. The MI
668     ! cell sizes are increased by 10% to provide a small buffer.
669           IF(USE_STL) THEN
670     
671              EXTENTS(1) = 6.0d0*MAX_DIA
672              EXTENTS(2) = WINDOW * 1.10d0
673              EXTENTS(3) = WINDOW * 1.10d0
674     
675              DO H=1,HMAX
676              DO W=1,WMAX
677     
678                 ORIGIN(1) = BC_X_w(BCV) - 3.0d0*MAX_DIA
679                 ORIGIN(2) = MESH_P(W) - WINDOW * 0.05d0
680                 ORIGIN(3) = MESH_Q(H) - WINDOW * 0.05d0
681     
682                 FACET_LP: DO LC=1, N_FACETS_DES
683     
684                    IF(BC_X_w(BCV) > maxval(VERTEX(:,1,LC))) CYCLE FACET_LP
685                    IF(BC_X_w(BCV) < minval(VERTEX(:,1,LC))) CYCLE FACET_LP
686     
687                    IF(BC_Y_s(BCV) > maxval(VERTEX(:,2,LC))) CYCLE FACET_LP
688                    IF(BC_Y_n(BCV) < minval(VERTEX(:,2,LC))) CYCLE FACET_LP
689     
690                    IF(BC_Z_b(BCV) > maxval(VERTEX(:,3,LC))) CYCLE FACET_LP
691                    IF(BC_Z_t(BCV) < minval(VERTEX(:,3,LC))) CYCLE FACET_LP
692     
693                    CALL TESTTRIANGLEAABB(VERTEX(:,:,LC), NORM_FACE(:,LC),  &
694                       ORIGIN(:), EXTENTS(:), SA_EXIST, SEP_AXIS, I, J, K)
695     
696                    IF(.NOT.SA_EXIST) THEN
697                       IF(NORM_FACE(2,LC) >= 0) THEN
698                          FULL_MAP(1:W,H) = 0
699                       ELSE
700                          FULL_MAP(W:WMAX,H) = 0
701                       ENDIF
702                       IF(NORM_FACE(3,LC) >= 0) THEN
703                          FULL_MAP(W,1:H) = 0
704                       ELSE
705                          FULL_MAP(W,H:HMAX) = 0
706                       ENDIF
707                    ENDIF
708                 ENDDO FACET_LP
709              ENDDO
710              ENDDO
711           ENDIF
712     
713     ! Add up the total number of available positions in the seeding map.
714           DO H=1,HMAX
715           DO W=1,WMAX
716              IF(FULL_MAP(W,H) /= 0) OCCUPANTS = OCCUPANTS + 1
717           ENDDO
718           ENDDO
719     
720     ! Sync the full map across all ranks.
721           CALL GLOBAL_ALL_SUM(OCCUPANTS)
722           CALL GLOBAL_ALL_SUM(FULL_MAP)
723     
724     ! Throw an error and exit if there are no occupants.
725           IF(OCCUPANTS == 0) THEN
726              WRITE(ERR_MSG, 1100) BCV_I
727              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
728           ENDIF
729     
730      1100 FORMAT('Error 1100: No un-cut fluid cells adjacent to DEM_MI ',  &
731              'staging area.',/'Unable to setup the discrete solids mass ', &
732              'inlet for BC:',I3)
733     
734     ! Store the number of occupants.
735           DEM_MI(BCV_I)%OCCUPANTS = OCCUPANTS
736     
737     ! Display the fill map if debugging
738           IF(dFlag .OR. showMAP) THEN
739              WRITE(*,"(2/,2x,'Displaying Fill Map:')")
740              DO H=HMAX,1,-1
741                 WRITE(*,"(2x,'H =',I3)",advance='no')H
742                 DO W=1,WMAX
743                    IF(FULL_MAP(W,H) == 0) then
744                       WRITE(*,"(' *')",advance='no')
745                    ELSE
746                       WRITE(*,"(' .')",advance='no')
747                    ENDIF
748                 ENDDO
749                 WRITE(*,*)' '
750              ENDDO
751           ENDIF
752     
753     ! Construct an array of integers randomly ordered.
754           if(dFLAG) write(*,"(2/,2x,'Building RAND_MAP:')")
755           allocate( RAND_MAP(OCCUPANTS) )
756           RAND_MAP = 0
757     
758     ! Only Rank 0 will randomize the layout and broadcast it to the other
759     ! ranks. This will ensure that all ranks see the same layout.
760           IF(myPE == 0) THEN
761              LL = 1
762              DO WHILE (RAND_MAP(OCCUPANTS) .EQ. 0)
763                 CALL RANDOM_NUMBER(TMP_DP)
764                 TMP_INT = CEILING(real(TMP_DP*dble(OCCUPANTS)))
765                 DO LC = 1, LL
766                   IF(TMP_INT .EQ. RAND_MAP(LC) )EXIT
767                   IF(LC .EQ. LL)THEN
768                      if(dFlag) WRITE(*,"(4x,'LC:',I6,' : ',I6)") LC, TMP_INT
769                      RAND_MAP(LC) = TMP_INT
770                      LL = LL + 1
771                   ENDIF
772                 ENDDO
773              ENDDO
774           ENDIF
775     
776           CALL GLOBAL_ALL_SUM(RAND_MAP)
777     
778     ! Initialize the vacancy pointer.
779           DEM_MI(BCV_I)%VACANCY = 1
780     
781     ! Allocate the mass inlet storage arrays.
782           allocate( DEM_MI(BCV_I)%W(OCCUPANTS) )
783           allocate( DEM_MI(BCV_I)%P(OCCUPANTS) )
784           allocate( DEM_MI(BCV_I)%H(OCCUPANTS) )
785           allocate( DEM_MI(BCV_I)%Q(OCCUPANTS) )
786           allocate( DEM_MI(BCV_I)%OWNER(OCCUPANTS) )
787     
788           if(dFlag) write(*,8010)
789     ! Store the MI layout in the randomized order.
790           LC = 0
791           DO H=1,HMAX
792           DO W=1,WMAX
793              IF(FULL_MAP(W,H) == 0) CYCLE
794              LC = LC + 1
795              LL = RAND_MAP(LC)
796              DEM_MI(BCV_I)%OWNER(LL) = FULL_MAP(W,H) - 1
797     
798              DEM_MI(BCV_I)%W(LL) = MESH_W(W)
799              DEM_MI(BCV_I)%H(LL) = MESH_H(H)
800     
801              DEM_MI(BCV_I)%P(LL) = MESH_P(W)
802              DEM_MI(BCV_I)%Q(LL) = MESH_Q(H)
803     
804              if(dFlag) write(*,8011) DEM_MI(BCV_I)%OWNER(LL), &
805                 DEM_MI(BCV_I)%W(LL), DEM_MI(BCV_I)%H(LL), DEM_MI(BCV_I)%L, &
806                 DEM_MI(BCV_I)%P(LL), DEM_MI(BCV_I)%Q(LL), DEM_MI(BCV_I)%OFFSET
807     
808           ENDDO
809           ENDDO
810     
811     
812      8010 FORMAT(2/,2x,'Storing DEM_MI data:',/4X,'OWNER',5X,'W',5X,'H',   &
813              5X,'L',7X,'P',12X,'Q',12X,'R')
814      8011 FORMAT(4x,I5,3(2X,I4),3(2x,g12.5))
815     
816     
817           if(dFlag) write(*,"(2/,2x,'Inlet area sizes:')")
818           if(dFlag) write(*,9000) 'mfix.dat: ', PLEN * QLEN
819           if(dFlag) write(*,9000) 'BC_AREA:  ', BC_AREA(BCV)
820           if(dFlag) write(*,9000) 'DEM_MI:   ', OCCUPANTS * (WINDOW**2)
821      9000 FORMAT(2x,A,g12.5)
822     
823     ! House keeping.
824           IF( allocated(MESH_H)) deallocate(MESH_H)
825           IF( allocated(MESH_W)) deallocate(MESH_W)
826           IF( allocated(MESH_P)) deallocate(MESH_P)
827           IF( allocated(MESH_Q)) deallocate(MESH_Q)
828     
829           IF( allocated(RAND_MAP)) deallocate(RAND_MAP)
830           IF( allocated(FULL_MAP)) deallocate(FULL_MAP)
831     
832           CALL FINL_ERR_MSG
833     
834           RETURN
835     
836      8005 FORMAT(2/,2x,'Building MESH_',A1,':',/4x,'Shift:',f8.4,/4x,      &
837              'MESH_',A1,'(0) = ',f8.4,/)
838     
839      8006 FORMAT(4x,'LC = ',I4,3x,A1,' =',I3,3x,A1,' =',f8.4)
840     
841           END SUBROUTINE LAYOUT_DEM_MI_EW
842     
843     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
844     !                                                                      !
845     !  Subroutine: LAYOUT_DEM_MI_TB                                        !
846     !                                                                      !
847     !  Purpose:  This routine determines the layout of the mass inlet as   !
848     !  either ordered or random based upon the inlet conditions.  This     !
849     !  routine also verifies that the specified inlet conditions for mass  !
850     !  or volumetric flow rates along with inlet size (length or area) and !
851     !  particle inlet velocity will work.  If not an error is flagged and  !
852     !  the program is exited.                                              !
853     !                                                                      !
854     !                                                                      !
855     !  Author: J.Musser                                   Date: 14-Aug-09  !
856     !                                                                      !
857     !  Comments:                                                           !
858     !                                                                      !
859     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
860           SUBROUTINE LAYOUT_DEM_MI_TB(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
861     
862     ! Global variables:
863     !---------------------------------------------------------------------//
864           use bc, only: BC_PLANE, BC_Z_b, BC_AREA
865           use bc, only: BC_X_w, BC_X_e
866           use bc, only: BC_Y_s, BC_Y_n
867     
868           use des_bc, only: DEM_MI
869     
870           use stl, only: N_FACETS_DES
871           use stl, only: VERTEX, NORM_FACE
872           use cutcell, only: USE_STL
873           use compar, only: myPE
874           use geometry, only: IMAX, JMAX, KMAX
875           use geometry, only: DX, DY, DZ
876           use geometry, only: XMIN, DO_K
877     
878           use funits, only: DMP_LOG
879     
880     ! Global parameters:
881     !---------------------------------------------------------------------//
882           use param1, only: ZERO, HALF, ONE
883     
884     ! Module procedures
885     !---------------------------------------------------------------------//
886           use des_bc, only: EXCLUDE_DEM_MI_CELL
887           use mpi_utility, only: GLOBAL_ALL_SUM
888           use des_stl_functions, only: TestTriangleAABB
889           use functions, only: IS_ON_myPE_OWNS
890     
891           use error_manager
892     
893           IMPLICIT NONE
894     
895     ! Dummy arguments
896     !---------------------------------------------------------------------//
897     ! passed arguments giving index information of the boundary
898           INTEGER, INTENT(IN) :: BCV
899           INTEGER, INTENT(IN) :: BCV_I
900     ! Max diameter of incoming particles at bc
901           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
902     ! Debug flas.
903           LOGICAL, INTENT(IN) :: setDBG, showMAP
904     
905     ! Local variables
906     !---------------------------------------------------------------------//
907     ! Loop counters.
908           INTEGER LL, LC
909     ! Indices for mapping to fluid grid.
910           INTEGER IJK, I, J, K
911     ! Local MESH indices
912           INTEGER H, W
913     ! Temp variable for double precision
914           DOUBLE PRECISION :: TMP_DP
915     ! Temporary variable for integers
916           INTEGER :: TMP_INT
917     ! Window indicies.
918           INTEGER, allocatable :: MESH_H(:)
919           INTEGER, allocatable :: MESH_W(:)
920     ! Window positions.
921           DOUBLE PRECISION, allocatable :: MESH_P(:)
922           DOUBLE PRECISION, allocatable :: MESH_Q(:)
923     ! Map of array index to MI cell
924           INTEGER, allocatable :: RAND_MAP(:)
925     ! Map of MI cells to owner processes
926           INTEGER, allocatable :: FULL_MAP(:,:)
927     ! max number of partitions along length of inlet
928           INTEGER :: WMAX, HMAX
929     ! the length of each side of the inlet boundary
930           DOUBLE PRECISION :: PLEN, QLEN
931     ! Number of occupied mesh cells
932           INTEGER :: OCCUPANTS
933     ! Offset and and window size.
934           DOUBLE PRECISION :: SHIFT, WINDOW
935     ! The origin and dimenion of MI cells. (STL intersection tests)
936           DOUBLE PRECISION :: ORIGIN(3), EXTENTS(3)
937     ! Separating axis test dummy variable
938           INTEGER :: SEP_AXIS
939     ! Indicates that a separating axis exists
940           LOGICAL :: SA_EXIST
941     ! Local Debug flag.
942           LOGICAL :: dFlag
943     
944     !......................................................................!
945     
946     ! Initialize the error manager.
947           CALL INIT_ERR_MSG('LAYOUT_DEM_MI_TB')
948     
949           dFlag = (DMP_LOG .AND. setDBG)
950           if(dFlag) write(*,"(2/,'Building TB DEM_MI: ',I3)") BCV_I
951     
952     ! Store the index that maps back to the user input.
953     
954           OCCUPANTS = 0
955     
956     ! Calculate number of partitions in first direction.
957           PLEN = BC_X_e(BCV) - BC_X_w(BCV)
958           WMAX = FLOOR(real(PLEN/MAX_DIA))
959           allocate( MESH_W(WMAX) )
960           allocate( MESH_P(0:WMAX) )
961     
962     ! Calculate number of partitions in the second direction.
963           QLEN = BC_Y_n(BCV) - BC_Y_s(BCV)
964           HMAX = FLOOR(real(QLEN/MAX_DIA))
965           allocate( MESH_H(HMAX) )
966           allocate( MESH_Q(0:HMAX) )
967     
968     ! Allocate the full map.
969           allocate( FULL_MAP(WMAX, HMAX))
970     
971     ! Set the value of the boundary condtion offset value used in the
972     ! placement of new particles.
973           CALL CALC_CELL_INTERSECT(ZERO, BC_Z_b(BCV), DZ, KMAX, K)
974           SHIFT = merge(-ONE, ONE, BC_PLANE(BCV) == 'T')
975           DEM_MI(BCV_I)%OFFSET = BC_Z_b(BCV) + MAX_DIA*SHIFT
976           DEM_MI(BCV_I)%L = K + int(SHIFT)
977           if(dFlag) write(*,"(2x,'Offset: ',3x,I4,3x,g12.5)") &
978              DEM_MI(BCV_I)%L, DEM_MI(BCV_I)%OFFSET
979     
980     
981     ! Dimension of grid cell for seeding particles; this may be larger than
982     ! than the particle diameter but not smaller:
983           DEM_MI(BCV_I)%WINDOW = MIN(PLEN/WMAX, QLEN/HMAX)
984           WINDOW = DEM_MI(BCV_I)%WINDOW
985           if(dFlag) write(*,"(2x,'Windows size: ',g12.5)") WINDOW
986     
987     ! Setup the first direction.
988           SHIFT = HALF*(PLEN - WMAX*WINDOW)
989           MESH_P(0) = BC_X_w(BCV) + SHIFT
990           if(dFlag) write(*,8005) 'P', SHIFT, 'P', MESH_P(0)
991           DO LC=1,WMAX
992              MESH_P(LC) = MESH_P(0) + dble(LC-1)*WINDOW
993              SHIFT = MESH_P(LC) + HALF*WINDOW
994              CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, MESH_W(LC))
995              IF(dFlag)WRITE(*,8006) LC, 'W', MESH_W(LC), 'P', MESH_P(LC)
996           ENDDO
997     
998     ! Setup the second direction.
999           SHIFT = HALF*(QLEN - HMAX*WINDOW)
1000           MESH_Q(0) = BC_Y_s(BCV) + SHIFT
1001           if(dFlag) write(*,8005) 'Q',SHIFT, 'Q',MESH_Q(0)
1002           DO LC=1,HMAX
1003              MESH_Q(LC) = MESH_Q(0) + dble(LC-1)*WINDOW
1004              SHIFT = MESH_Q(LC) + HALF*WINDOW
1005              CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, MESH_H(LC))
1006              IF(dFlag)WRITE(*,8006) LC, 'H', MESH_H(LC), 'Q', MESH_Q(LC)
1007           ENDDO
1008     
1009     ! Get the Jth index of the fluid cell
1010           CALL CALC_CELL_INTERSECT(ZERO, BC_Z_b(BCV), DZ, KMAX, K)
1011     
1012     ! If the computationsl cell adjacent to the DEM_MI mesh cell is a
1013     ! fluid cell and has not been cut, store the ID of the cell owner.
1014           DO H=1,HMAX
1015           DO W=1,WMAX
1016     
1017              I = MESH_W(W)
1018              J = MESH_H(H)
1019              FULL_MAP(W,H) = 0
1020     
1021              IF(.NOT.IS_ON_myPE_owns(I,J,K)) CYCLE
1022     
1023              CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
1024              CALL CALC_CELL_INTERSECT(ZERO, MESH_Q(H), DY, JMAX, J)
1025              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1026     
1027              SHIFT = MESH_P(W)+WINDOW
1028              CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, I)
1029              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1030     
1031              SHIFT = MESH_Q(H)+WINDOW
1032              CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, J)
1033              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1034     
1035              CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
1036              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1037     
1038              FULL_MAP(W,H) = myPE+1
1039           ENDDO
1040           ENDDO
1041     
1042     ! For complex boundaries defined by STLs, exclude any mass inflow cells
1043     ! that intersect the boundary and all cells opposite the normal. The MI
1044     ! cell sizes are increased by 10% to provide a small buffer.
1045           IF(USE_STL) THEN
1046     
1047              EXTENTS(3) = 6.0d0*MAX_DIA
1048              EXTENTS(1) = WINDOW * 1.10d0
1049              EXTENTS(2) = WINDOW * 1.10d0
1050     
1051              DO H=1,HMAX
1052              DO W=1,WMAX
1053     
1054                 ORIGIN(3) = BC_Z_b(BCV) - 3.0d0*MAX_DIA
1055                 ORIGIN(1) = MESH_P(W) - WINDOW * 0.05d0
1056                 ORIGIN(2) = MESH_Q(H) - WINDOW * 0.05d0
1057     
1058                 FACET_LP: DO LC=1, N_FACETS_DES
1059     
1060                    IF(BC_Z_b(BCV) > maxval(VERTEX(:,3,LC))) CYCLE FACET_LP
1061                    IF(BC_Z_b(BCV) < minval(VERTEX(:,3,LC))) CYCLE FACET_LP
1062     
1063                    IF(BC_X_w(BCV) > maxval(VERTEX(:,1,LC))) CYCLE FACET_LP
1064                    IF(BC_X_e(BCV) < minval(VERTEX(:,1,LC))) CYCLE FACET_LP
1065     
1066                    IF(BC_Y_s(BCV) > maxval(VERTEX(:,2,LC))) CYCLE FACET_LP
1067                    IF(BC_Y_n(BCV) < minval(VERTEX(:,2,LC))) CYCLE FACET_LP
1068     
1069                    CALL TESTTRIANGLEAABB(VERTEX(:,:,LC), NORM_FACE(:,LC),  &
1070                       ORIGIN(:), EXTENTS(:), SA_EXIST, SEP_AXIS, I, J, K)
1071     
1072                    IF(.NOT.SA_EXIST) THEN
1073                       IF(NORM_FACE(1,LC) >= 0) THEN
1074                          FULL_MAP(1:W,H) = 0
1075                       ELSE
1076                          FULL_MAP(W:WMAX,H) = 0
1077                       ENDIF
1078                       IF(NORM_FACE(2,LC) >= 0) THEN
1079                          FULL_MAP(W,1:H) = 0
1080                       ELSE
1081                          FULL_MAP(W,H:HMAX) = 0
1082                       ENDIF
1083                    ENDIF
1084                 ENDDO FACET_LP
1085              ENDDO
1086              ENDDO
1087           ENDIF
1088     
1089     ! Add up the total number of available positions in the seeding map.
1090           DO H=1,HMAX
1091           DO W=1,WMAX
1092              IF(FULL_MAP(W,H) /= 0) OCCUPANTS = OCCUPANTS + 1
1093           ENDDO
1094           ENDDO
1095     
1096     ! Sync the full map across all ranks.
1097           CALL GLOBAL_ALL_SUM(OCCUPANTS)
1098           CALL GLOBAL_ALL_SUM(FULL_MAP)
1099     
1100     ! Throw an error and exit if there are no occupants.
1101           IF(OCCUPANTS == 0) THEN
1102              WRITE(ERR_MSG, 1100) BCV_I
1103              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
1104           ENDIF
1105     
1106      1100 FORMAT('Error 1100: No un-cut fluid cells adjacent to DEM_MI ',  &
1107              'staging area.',/'Unable to setup the discrete solids mass ', &
1108              'inlet for BC:',I3)
1109     
1110     ! Store the number of occupants.
1111           DEM_MI(BCV_I)%OCCUPANTS = OCCUPANTS
1112     
1113     ! Display the fill map if debugging
1114           IF(dFlag .OR. showMAP) THEN
1115              WRITE(*,"(2/,2x,'Displaying Fill Map:')")
1116              DO H=HMAX,1,-1
1117                 WRITE(*,"(2x,'H =',I3)",advance='no')H
1118                 DO W=1,WMAX
1119                    IF(FULL_MAP(W,H) == 0) then
1120                       WRITE(*,"(' *')",advance='no')
1121                    ELSE
1122                       WRITE(*,"(' .')",advance='no')
1123                    ENDIF
1124                 ENDDO
1125                 WRITE(*,*)' '
1126              ENDDO
1127           ENDIF
1128     
1129     ! Construct an array of integers randomly ordered.
1130           if(dFLAG) write(*,"(2/,2x,'Building RAND_MAP:')")
1131           allocate( RAND_MAP(OCCUPANTS) )
1132           RAND_MAP = 0
1133     
1134     ! Only Rank 0 will randomize the layout and broadcast it to the other
1135     ! ranks. This will ensure that all ranks see the same layout.
1136           IF(myPE == 0) THEN
1137              LL = 1
1138              DO WHILE (RAND_MAP(OCCUPANTS) .EQ. 0)
1139                 CALL RANDOM_NUMBER(TMP_DP)
1140                 TMP_INT = CEILING(real(TMP_DP*dble(OCCUPANTS)))
1141                 DO LC = 1, LL
1142                   IF(TMP_INT .EQ. RAND_MAP(LC) )EXIT
1143                   IF(LC .EQ. LL)THEN
1144                      if(dFlag) WRITE(*,"(4x,'LC:',I3,' : ',I3)") LC, TMP_INT
1145                      RAND_MAP(LC) = TMP_INT
1146                      LL = LL + 1
1147                   ENDIF
1148                 ENDDO
1149              ENDDO
1150           ENDIF
1151     
1152           CALL GLOBAL_ALL_SUM(RAND_MAP)
1153     
1154     ! Initialize the vacancy pointer.
1155           DEM_MI(BCV_I)%VACANCY = 1
1156     
1157     ! Allocate the mass inlet storage arrays.
1158           allocate( DEM_MI(BCV_I)%W(OCCUPANTS) )
1159           allocate( DEM_MI(BCV_I)%P(OCCUPANTS) )
1160           allocate( DEM_MI(BCV_I)%H(OCCUPANTS) )
1161           allocate( DEM_MI(BCV_I)%Q(OCCUPANTS) )
1162           allocate( DEM_MI(BCV_I)%OWNER(OCCUPANTS) )
1163     
1164           if(dFlag) write(*,8010)
1165     ! Store the MI layout in the randomized order.
1166           LC = 0
1167           DO H=1,HMAX
1168           DO W=1,WMAX
1169              IF(FULL_MAP(W,H) == 0) CYCLE
1170              LC = LC + 1
1171              LL = RAND_MAP(LC)
1172              DEM_MI(BCV_I)%OWNER(LL) = FULL_MAP(W,H) - 1
1173     
1174              DEM_MI(BCV_I)%W(LL) = MESH_W(W)
1175              DEM_MI(BCV_I)%H(LL) = MESH_H(H)
1176     
1177              DEM_MI(BCV_I)%P(LL) = MESH_P(W)
1178              DEM_MI(BCV_I)%Q(LL) = MESH_Q(H)
1179     
1180              if(dFlag) write(*,8011) DEM_MI(BCV_I)%OWNER(LL), &
1181                 DEM_MI(BCV_I)%W(LL), DEM_MI(BCV_I)%H(LL), DEM_MI(BCV_I)%L, &
1182                 DEM_MI(BCV_I)%P(LL), DEM_MI(BCV_I)%Q(LL), DEM_MI(BCV_I)%OFFSET
1183     
1184           ENDDO
1185           ENDDO
1186     
1187     
1188      8010 FORMAT(2/,2x,'Storing DEM_MI data:',/4X,'OWNER',5X,'W',5X,'H',   &
1189              5X,'L',7X,'P',12X,'Q',12X,'R')
1190      8011 FORMAT(4x,I5,3(2X,I4),3(2x,g12.5))
1191     
1192     
1193           if(dFlag) write(*,"(2/,2x,'Inlet area sizes:')")
1194           if(dFlag) write(*,9000) 'mfix.dat: ', PLEN * QLEN
1195           if(dFlag) write(*,9000) 'BC_AREA:  ', BC_AREA(BCV)
1196           if(dFlag) write(*,9000) 'DEM_MI:   ', OCCUPANTS * (WINDOW**2)
1197      9000 FORMAT(2x,A,g12.5)
1198     
1199     ! House keeping.
1200           IF( allocated(MESH_H)) deallocate(MESH_H)
1201           IF( allocated(MESH_W)) deallocate(MESH_W)
1202           IF( allocated(MESH_P)) deallocate(MESH_P)
1203           IF( allocated(MESH_Q)) deallocate(MESH_Q)
1204     
1205           IF( allocated(RAND_MAP)) deallocate(RAND_MAP)
1206           IF( allocated(FULL_MAP)) deallocate(FULL_MAP)
1207     
1208           CALL FINL_ERR_MSG
1209     
1210           RETURN
1211     
1212      8005 FORMAT(2/,2x,'Building MESH_',A1,':',/4x,'Shift:',f8.4,/4x,      &
1213              'MESH_',A1,'(0) = ',f8.4,/)
1214     
1215      8006 FORMAT(4x,'LC = ',I4,3x,A1,' =',I3,3x,A1,' =',f8.4)
1216     
1217           END SUBROUTINE LAYOUT_DEM_MI_TB
1218     
1219     
1220     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1221     !                                                                      !
1222     !  Subroutine: SET_DEM_MI_OWNER                                        !
1223     !  Author: J.Musser                                   Date: 08-Aug-14  !
1224     !                                                                      !
1225     !  Purpose: Set the owner of the background mesh used for seeding new  !
1226     !  particles into the domain. The mesh is stored and read from the RES !
1227     !  file, however, the owers need to be set per-run in case the DMP     !
1228     !  partitions changed.                                                 !
1229     !                                                                      !
1230     !  Comments:                                                           !
1231     !                                                                      !
1232     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1233           SUBROUTINE SET_DEM_MI_OWNER(BCV, BCV_I)
1234     
1235           use bc, only: BC_PLANE
1236           use des_bc, only: DEM_MI
1237     
1238           use compar
1239           use geometry
1240           use indices
1241     
1242           use funits, only: DMP_LOG
1243     
1244     ! Module procedures
1245     !---------------------------------------------------------------------//
1246           use mpi_utility, only: GLOBAL_ALL_SUM
1247           use error_manager
1248           use functions
1249     
1250           IMPLICIT NONE
1251     !-----------------------------------------------
1252     ! Dummy arguments
1253     !-----------------------------------------------
1254     ! passed arguments giving index information of the boundary
1255           INTEGER, INTENT(IN) :: BCV
1256           INTEGER, INTENT(IN) :: BCV_I
1257     
1258     ! Generic loop counter
1259           INTEGER :: LC1, OCCUPANTS
1260           INTEGER :: I, J, K
1261     
1262           OCCUPANTS = DEM_MI(BCV_I)%OCCUPANTS
1263           allocate(DEM_MI(BCV_I)%OWNER(OCCUPANTS))
1264     
1265           DEM_MI(BCV_I)%OWNER = 0
1266     
1267           SELECT CASE (BC_PLANE(BCV))
1268           CASE('N','S')
1269     
1270              J = DEM_MI(BCV_I)%L
1271              DO LC1=1,OCCUPANTS
1272                 I = DEM_MI(BCV_I)%W(LC1)
1273                 K = DEM_MI(BCV_I)%H(LC1)
1274                 IF(IS_ON_myPE_owns(I,J,K)) &
1275                    DEM_MI(BCV_I)%OWNER(LC1) = myPE
1276              ENDDO
1277     
1278           CASE('E','W')
1279     
1280              I = DEM_MI(BCV_I)%L
1281              DO LC1=1,OCCUPANTS
1282                 J = DEM_MI(BCV_I)%W(LC1)
1283                 K = DEM_MI(BCV_I)%H(LC1)
1284                 IF(IS_ON_myPE_owns(I,J,K)) &
1285                    DEM_MI(BCV_I)%OWNER(LC1) = myPE
1286              ENDDO
1287     
1288           CASE('T','B')
1289     
1290              K = DEM_MI(BCV_I)%L
1291              DO LC1=1,OCCUPANTS
1292                 I = DEM_MI(BCV_I)%W(LC1)
1293                 J = DEM_MI(BCV_I)%H(LC1)
1294                 IF(IS_ON_myPE_owns(I,J,K)) &
1295                    DEM_MI(BCV_I)%OWNER(LC1) = myPE
1296              ENDDO
1297     
1298           END SELECT
1299     
1300           CALL GLOBAL_ALL_SUM(DEM_MI(BCV_I)%OWNER(:))
1301     
1302           RETURN
1303           END SUBROUTINE SET_DEM_MI_OWNER
1304