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

1           MODULE vtp
2     
3           use mpi_utility
4           use cdist
5     
6           use desmpi
7           use mpi_comm_des
8           use error_manager
9     
10           IMPLICIT NONE
11     
12           INTEGER, PRIVATE :: GLOBAL_CNT
13           INTEGER, PRIVATE :: LOCAL_CNT
14     
15           INTEGER :: DES_UNIT = 2000
16     
17     ! file unit for ParaView *.pvd data
18           INTEGER, PARAMETER :: PVD_UNIT = 2050
19     
20     ! formatted file name
21           CHARACTER(LEN=511) :: FNAME_VTP
22     
23           INTERFACE VTP_WRITE_DATA
24              MODULE PROCEDURE VTP_WRITE_DP1
25              MODULE PROCEDURE VTP_WRITE_DP2
26              MODULE PROCEDURE VTP_WRITE_I1
27           END INTERFACE
28     
29           CONTAINS
30     
31     !``````````````````````````````````````````````````````````````````````!
32     ! Subroutine: VTP_WRITE_DP1                                            !
33     !                                                                      !
34     ! Purpose: Collect and write 1D double percision arrays to the VTP     !
35     ! file. This routine is designed to collect the data for parallel and  !
36     ! serial runs. This routine also manages the distribted IO case.       !
37     !``````````````````````````````````````````````````````````````````````!
38           SUBROUTINE VTP_WRITE_DP1(NAME, DATA)
39     
40           CHARACTER(len=*), INTENT(in) :: NAME
41           DOUBLE PRECISION, INTENT(in) :: DATA(:)
42     
43           INTEGER :: LC, PC
44     
45           IF(bDist_IO) THEN
46     
47              WRITE(DES_UNIT,1000) NAME
48     
49              PC = 1
50              DO LC = 1, MAX_PIP
51                 IF(PC > PIP) EXIT
52                 IF(IS_NONEXISTENT(LC)) CYCLE
53                 PC = PC+1
54                 IF(IS_GHOST(LC) .OR. IS_ENTERING_GHOST(LC) .OR. IS_EXITING_GHOST(LC)) CYCLE
55                 WRITE(DES_UNIT, 1001,ADVANCE="NO") real(DATA(LC))
56              ENDDO
57              WRITE(DES_UNIT,1002)
58     
59           ELSE
60     
61              allocate (dProcBuf(LOCAL_CNT) )
62              allocate (dRootBuf(GLOBAL_CNT))
63     
64              CALL DES_GATHER(DATA)
65     
66              IF(myPE == PE_IO) THEN
67                 WRITE(DES_UNIT,1000) NAME
68                 DO LC=1, GLOBAL_CNT
69                    WRITE(DES_UNIT,1001,ADVANCE="NO") real(drootbuf(LC))
70                 ENDDO
71                 WRITE(DES_UNIT,1002)
72              ENDIF
73     
74              deallocate(dProcBuf, dRootBuf)
75     
76           ENDIF
77     
78      1000 FORMAT('<DataArray type="Float32" Name="',A,'" format="ascii">')
79      1001 FORMAT(ES14.6,1X)
80      1002 FORMAT('</DataArray>')
81     
82           END SUBROUTINE VTP_WRITE_DP1
83     
84     !``````````````````````````````````````````````````````````````````````!
85     ! Subroutine: VTP_WRITE_DP2                                            !
86     !                                                                      !
87     ! Purpose: Collect and write 2D double percision arrays to the VTP     !
88     ! file. This routine is designed to collect the data for parallel and  !
89     ! serial runs. This routine also manages the distribted IO case.       !
90     !``````````````````````````````````````````````````````````````````````!
91           SUBROUTINE VTP_WRITE_DP2(NAME, DATA)
92     
93           CHARACTER(len=*), INTENT(in) :: NAME
94           DOUBLE PRECISION, INTENT(in) :: DATA(:,:)
95     
96           DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)
97     
98           CHARACTER(len=16) :: NOC
99           INTEGER :: LB, UB
100           INTEGER :: PC, LC1, LC2
101     
102           LB = LBOUND(DATA,2)
103           UB = UBOUND(DATA,2)
104           NOC=''; WRITE(NOC,*) (UB-LB)+1
105     
106           IF(bDist_IO) THEN
107     
108              WRITE(DES_UNIT,1000) NAME, trim(adjustl(NOC))
109     
110              PC = 1
111              DO LC1 = 1, MAX_PIP
112                 IF(PC > PIP) EXIT
113                 IF(IS_NONEXISTENT(LC1)) CYCLE
114                 PC = PC+1
115                 IF(IS_GHOST(LC1) .OR. IS_ENTERING_GHOST(LC1) .OR. IS_EXITING_GHOST(LC1)) CYCLE
116                 DO LC2=LB, UB
117                    WRITE(DES_UNIT,1001,ADVANCE="NO") real(DATA(LC1,LC2))
118                 ENDDO
119              ENDDO
120              WRITE(DES_UNIT,1002)
121     
122           ELSE
123     
124              allocate (dProcBuf(LOCAL_CNT) )
125              allocate (dRootBuf(GLOBAL_CNT))
126              allocate (ltemp_array((UB-LB)+1,GLOBAL_CNT))
127     
128              DO LC1 = LB, UB
129                 CALL DES_GATHER(DATA(:,LC1))
130                 ltemp_array(LC1,:) = drootbuf(:)
131              ENDDO
132     
133              IF(myPE == PE_IO) THEN
134                 WRITE(DES_UNIT,1000) NAME, trim(adjustl(NOC))
135                 DO LC1=1, GLOBAL_CNT
136                    DO LC2=LB, UB
137                       WRITE(DES_UNIT,1001,ADVANCE="NO") &
138                          real(ltemp_array(LC2,LC1))
139                    ENDDO
140                 ENDDO
141                 WRITE(DES_UNIT,1002)
142              ENDIF
143     
144              deallocate (dProcBuf, dRootBuf, ltemp_array)
145     
146           ENDIF
147     
148     
149      1000 FORMAT('<DataArray type="Float32" Name="',A,'" NumberOf',        &
150              'Components="',A,'" format="ascii">')
151      1001 FORMAT(ES14.6,1X)
152      1002 FORMAT('</DataArray>')
153     
154           END SUBROUTINE VTP_WRITE_DP2
155     
156     
157     
158     !``````````````````````````````````````````````````````````````````````!
159     ! Subroutine: VTP_WRITE_I1                                             !
160     !                                                                      !
161     ! Purpose: Collect and write 1D integer arrays to the VTP file. This   !
162     ! routine is designed to collect the data for parallel and serial      !
163     ! runs. This routine also manages the distribted IO case.              !
164     !``````````````````````````````````````````````````````````````````````!
165           SUBROUTINE VTP_WRITE_I1(NAME, DATA)
166     
167           CHARACTER(len=*), INTENT(in) :: NAME
168           INTEGER, INTENT(in) :: DATA(:)
169     
170           INTEGER :: LC, PC
171     
172           IF(bDist_IO) THEN
173     
174              WRITE(DES_UNIT,1000) NAME
175     
176              PC = 1
177              DO LC = 1, MAX_PIP
178                 IF(PC > PIP) EXIT
179                 IF(IS_NONEXISTENT(LC)) CYCLE
180                 PC = PC+1
181                 IF(IS_GHOST(LC) .OR. IS_ENTERING_GHOST(LC) .OR. IS_EXITING_GHOST(LC)) CYCLE
182                 WRITE(DES_UNIT, 1001,ADVANCE="NO") DATA(LC)
183              ENDDO
184              WRITE(DES_UNIT,1002)
185     
186           ELSE
187     
188              allocate (iProcBuf(LOCAL_CNT) )
189              allocate (iRootBuf(GLOBAL_CNT))
190     
191              CALL DES_GATHER(DATA)
192     
193              IF(myPE == PE_IO) THEN
194                 WRITE(DES_UNIT,1000) NAME
195                 DO LC=1, GLOBAL_CNT
196                    WRITE(DES_UNIT,1001,ADVANCE="NO") irootbuf(LC)
197                 ENDDO
198                 WRITE(DES_UNIT,1002)
199              ENDIF
200     
201              deallocate(iProcBuf, iRootBuf)
202     
203           ENDIF
204     
205      1000 FORMAT('<DataArray type="Float32" Name="',A,'" format="ascii">')
206      1001 FORMAT(I10,1X)
207      1002 FORMAT('</DataArray>')
208     
209           END SUBROUTINE VTP_WRITE_I1
210     
211     
212     !``````````````````````````````````````````````````````````````````````!
213     ! Subroutine: VTP_WRITE_ELEMENT                                        !
214     !                                                                      !
215     ! Purpose: Write a string to the VTP file. It masks the need to check  !
216     ! the logical before flushing.                                         !
217     !``````````````````````````````````````````````````````````````````````!
218           SUBROUTINE VTP_WRITE_ELEMENT(ELEMENT)
219     
220           CHARACTER(len=*), INTENT(in) :: ELEMENT
221     
222           IF(bDist_IO .OR. myPE == PE_IO) &
223              WRITE(DES_UNIT,"(A)") ELEMENT
224     
225           RETURN
226           END SUBROUTINE VTP_WRITE_ELEMENT
227     
228     
229     
230     !``````````````````````````````````````````````````````````````````````!
231     ! Subroutine: VTP_OPEN_FILE                                            !
232     !                                                                      !
233     ! Purpose: This routine opens the VTP file and calcualtes the offsets  !
234     ! for dmp data collection.                                             !
235     !``````````````````````````````````````````````````````````````````````!
236           SUBROUTINE VTP_OPEN_FILE(NoPc)
237     
238     ! Modules
239     !-----------------------------------------------
240           use discretelement, only: VTP_DIR
241           use run, only: RUN_TYPE, RUN_NAME
242     
243           IMPLICIT NONE
244     
245           CHARACTER(len=*) :: NoPc
246     
247           INTEGER :: NumberOfPoints
248     
249     ! Variables related to gather
250           integer lgathercnts(0:numpes-1), lproc
251     
252     ! check whether an error occurs in opening a file
253           INTEGER :: IOS
254     ! Integer error flag.
255           INTEGER :: IER
256     
257     ! logical used for testing is the data file already exists
258           LOGICAL :: EXISTS_VTP
259     ! status of the vtp file to be written
260           CHARACTER(LEN=8) :: STATUS_VTP
261     
262           IF(TRIM(VTP_DIR)/='.') CALL CREATE_DIR(trim(VTP_DIR))
263     
264     ! Initial the global count.
265           GLOBAL_CNT = 10
266     ! Calculate the number of 'real' particles on the local process.
267           LOCAL_CNT = PIP - iGHOST_CNT
268     
269     ! Distributed IO
270           IF(bDIST_IO) THEN
271              NumberOfPoints = LOCAL_CNT
272              WRITE(NoPc,"(I10.10)") NumberOfPoints
273     
274              IF(TRIM(VTP_DIR)/='.') THEN
275                 WRITE(fname_vtp,'(A,"/",A,"_DES",I4.4,"_",I5.5,".vtp")') &
276                    trim(VTP_DIR), trim(run_name), vtp_findex, mype
277              ELSE
278                 WRITE(fname_vtp,'(A,"_DES",I4.4,"_",I5.5,".vtp")') &
279                    trim(run_name), vtp_findex, mype
280              ENDIF
281     
282     ! Serial IO
283           ELSE
284     
285     ! Calculate the total number of particles system-wide.
286              call global_sum(LOCAL_CNT, GLOBAL_CNT)
287              NumberOfPoints = GLOBAL_CNT
288              WRITE(NoPc,"(I10.10)") NumberOfPoints
289     
290     ! Set the send count from the local process.
291              igath_sendcnt = LOCAL_CNT
292     
293     ! Collect the number of particles on each rank.all ranks.
294              lgathercnts = 0
295              lgathercnts(myPE) = LOCAL_CNT
296              call global_sum(lgathercnts,igathercnts)
297     
298     ! Calculate the rank displacements.
299              idispls(0) = 0
300              DO lPROC = 1,NUMPEs-1
301                 idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
302              ENDDO
303     
304     ! set the file name and unit number and open file
305              IF(TRIM(VTP_DIR)/='.') THEN
306                 WRITE(fname_vtp,'(A,"/",A,"_DES_",I5.5,".vtp")') &
307                    trim(VTP_DIR),trim(run_name), vtp_findex
308              ELSE
309                 WRITE(fname_vtp,'(A,"_DES_",I5.5,".vtp")') &
310                    trim(run_name), vtp_findex
311              ENDIF
312           ENDIF
313     
314           IER = 0
315           IF(bDIST_IO .OR. myPE == PE_IO) THEN
316     
317     ! The file should be new but could exist due to restarting.
318              STATUS_VTP = 'NEW'
319     ! Check to see if the file already exists.
320              INQUIRE(FILE=FNAME_VTP,EXIST=EXISTS_VTP)
321     ! The given file should not exist if the run type is NEW.
322              IF(EXISTS_VTP)THEN
323     ! The VTP should never exist for a NEW run.
324                 IF(RUN_TYPE == 'NEW')THEN
325                    IER = 1
326     ! The file may exist during a RESTART.
327                 ELSE
328                    STATUS_VTP = 'REPLACE'
329                 ENDIF
330              ENDIF
331     
332     ! Open the file and record any erros.
333              IF(IER == 0) THEN
334                 OPEN(CONVERT='BIG_ENDIAN',UNIT=DES_UNIT, FILE=FNAME_VTP,   &
335                    STATUS=STATUS_VTP, IOSTAT=IOS)
336                 IF(IOS /= 0) IER = 2
337              ENDIF
338           ENDIF
339     
340           CALL GLOBAL_ALL_MAX(IER)
341     
342           IF(IER /= 0) THEN
343              CALL INIT_ERR_MSG("VTP_MOD --> OPEN_VTP")
344              WRITE(ERR_MSG, 1100) IER
345              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
346           ENDIF
347     
348      1100 FORMAT('Error 1100: Unable to open VTP file. This could be ',    &
349              'caused by a VTP',/'file with the same file name already ',   &
350              'existing. or an error code',/' returned by the OPEN ',       &
351              'function.'/'Error code: ',I2,4x,'Aborting.')
352     
353     
354           END SUBROUTINE VTP_OPEN_FILE
355     
356     
357     
358     !......................................................................!
359     ! SUBROUTINE: VTP_CLOSE_FILE                                           !
360     !                                                                      !
361     ! Purpose: This routine closes the vtp file.                           !
362     !``````````````````````````````````````````````````````````````````````!
363           SUBROUTINE VTP_CLOSE_FILE
364     
365     
366           VTP_FINDEX=VTP_FINDEX+1
367     
368           IF(bDist_io .OR. (myPE .eq.pe_IO)) CLOSE(des_unit)
369     
370     
371           END SUBROUTINE VTP_CLOSE_FILE
372     
373     
374     !......................................................................!
375     ! SUBROUTINE: ADD_VTP_TO_PVD                                           !
376     !                                                                      !
377     ! Purpose: This routine opens the pvd file.                            !
378     !``````````````````````````````````````````````````````````````````````!
379           SUBROUTINE ADD_VTP_TO_PVD
380     
381           use discretelement, only: VTP_DIR
382           use run, only: RUN_TYPE, RUN_NAME
383     
384     !-----------------------------------------------
385     ! Local Variables
386     !-----------------------------------------------
387     ! Index position of desired character
388           INTEGER IDX_f, IDX_b
389     ! logical used for testing is the data file already exists
390           LOGICAL :: EXISTS_PVD
391     ! Generic input limited to 256 characters
392           CHARACTER(LEN=256) INPUT
393     
394     ! formatted file name
395           CHARACTER(LEN=64) :: FNAME_PVD = ''
396     ! formatted time
397           CHARACTER(LEN=64) :: cTIME = ''
398     
399           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
400     
401     ! IO Status flag
402           INTEGER :: IOS
403     
404     ! Variables related to gather
405           integer :: IER
406     
407     !-----------------------------------------------
408     
409           CALL INIT_ERR_MSG('VTP_MOD --> ADD_VTP_TO_PVD')
410     
411     ! Initialize the error flag.
412           IER = 0
413     
414     ! Obtain the file name and open the pvd file
415           FNAME_PVD = TRIM(RUN_NAME)//'_DES.pvd'
416     
417     ! The PVD file is only written by PE_IO with serial IO.
418           IF(myPE == PE_IO .AND. .NOT.bDist_IO) THEN
419     
420     ! Check to see if the file already exists.
421              INQUIRE(FILE=FNAME_PVD,EXIST=EXISTS_PVD)
422     
423              IF(FIRST_PASS) THEN
424     
425     ! Open the "NEW" file and write the necessary header information.
426                 IF(RUN_TYPE /= 'RESTART_1')THEN
427     
428     ! The file exists but first_pass is also true so most likely an existing
429     ! file from an earlier/other run is present in the directory. Exit to
430     ! prevent accidently overwriting the existing file.
431                    IF(EXISTS_PVD) THEN
432                       IER = 1
433                    ELSE
434                       OPEN(UNIT=PVD_UNIT,FILE=FNAME_PVD,STATUS='NEW')
435                       WRITE(PVD_UNIT,"(A)")'<?xml version="1.0"?>'
436                       WRITE(PVD_UNIT,"(A)")'<VTKFile type="Collection" &
437                          &version="0.1" byte_order="LittleEndian">'
438                       WRITE(PVD_UNIT,"(3X,'<Collection>')")
439                    ENDIF
440     
441     ! This is the first pass of a restart run. Extra care is needed to make
442     ! sure that the pvd file is ready to accept new data.
443                 ELSE ! a restart run
444                    IF(EXISTS_PVD) THEN
445     ! Open the file at the beginning.
446                       OPEN(UNIT=PVD_UNIT,FILE=FNAME_PVD,&
447                          POSITION="REWIND",STATUS='OLD',IOSTAT=IOS)
448                       IF(IOS /= 0) IER = 2
449                    ELSE ! a pvd file does not exist
450                       IER = 3
451                    ENDIF
452     
453                    IF(IER == 0) THEN
454     ! Loop over the entries in the PVD file, looking for a match to the
455     ! file that is being written. If no match is found, the data will be
456     ! appended to the end of the pvd file, otherwise, the old data will
457     ! be over-written.
458                       DO
459     ! Read in the entires of the PVD file.
460                          READ(PVD_UNIT,"(A)",IOSTAT=IOS)INPUT
461                          IF(IOS > 0) THEN
462                             IER = 4
463                             EXIT
464                          ELSEIF(IOS<0)THEN
465     ! The end of the pvd file has been reached without finding an entry
466     ! matching the current record. Exit the loop.
467                             BACKSPACE(PVD_UNIT)
468                             EXIT
469                          ENDIF
470     ! Find the first instances of file=" and "/> in the read data.
471                          IDX_f = INDEX(INPUT,'file="')
472                          IDX_b = INDEX(INPUT,'"/>')
473     ! Skip rows that do not contain file data
474                          IF(IDX_f == 0 .AND. IDX_b == 0) CYCLE
475     ! Truncate the file name from the read data
476                          WRITE (INPUT,"(A)") INPUT(IDX_f+6:IDX_b-1)
477     ! If the file name matches the current VTP record, break the loop to
478     ! over-write this record.
479                          IF(TRIM(FNAME_VTP) == TRIM(INPUT)) THEN
480                             BACKSPACE(PVD_UNIT)
481                             EXIT
482                          ENDIF
483                       ENDDO
484                    ENDIF ! No errors
485                 ENDIF ! run_type new or restart
486     
487              ELSE ! not FIRST_PASS
488                 OPEN(UNIT=PVD_UNIT,FILE=FNAME_PVD,&
489                    POSITION="APPEND",STATUS='OLD',IOSTAT=IOS)
490                 IF (IOS /= 0) IER = 2
491              ENDIF
492     
493           ENDIF ! if myPE == PE_IO and not distributed IO
494     
495     
496           CAlL GLOBAL_ALL_SUM(IER)
497           IF(IER /= 0) THEN
498              SELECT CASE(IER)
499              CASE(1); WRITE(ERR_MSG,1101) trim(FNAME_PVD)
500              CASE(2); WRITE(ERR_MSG,1102) trim(FNAME_PVD)
501              CASE(3); WRITE(ERR_MSG,1103) trim(FNAME_PVD)
502              CASE(4); WRITE(ERR_MSG,1104) trim(FNAME_PVD)
503              CASE DEFAULT; WRITE(ERR_MSG,1105) trim(FNAME_PVD)
504              END SELECT
505              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
506           ENDIF
507     
508      1101 FORMAT('Error 1101: A PVD file was detected in the run ',        &
509              'directory which should',/'not exist for a NEW run.',/        &
510              'File: ',A)
511     
512      1102 FORMAT('Error 1102: Fatal error status returned while OPENING ', &
513              'PVD file.',/'File: ', A)
514     
515      1103 FORMAT('Error 1103: PVD file MISSING from run directory.',/      &
516              'File: ',A)
517     
518      1104 FORMAT('Error 1104: Fatal error status returned while READING ', &
519              'PVD file.',/'File: ', A)
520     
521      1105 FORMAT('Error 1105:: Fatal unclassified error when processing ', &
522              'PVD file.',/'File: ', A)
523     
524     
525     ! If there were no errors, updated the file.
526           IF(myPE == PE_IO .AND. .NOT.bDist_IO) THEN
527     
528     ! Remove the last two lines written so that additional data can be added
529              IF(.NOT.FIRST_PASS) THEN
530                 BACKSPACE(PVD_UNIT)
531                 BACKSPACE(PVD_UNIT)
532              ENDIF
533     
534              WRITE(cTIME,"(F12.6)") S_TIME
535     ! Write the data to the file
536              WRITE(PVD_UNIT,"(6X,A,A,A,A,A,A,A)")&
537              '<DataSet timestep="',trim(adjustl(cTIME)),'" ',&
538              'group="" part="0" ',& ! necessary file data
539              'file="',TRIM(FNAME_VTP),'"/>' ! file name of vtp
540     
541     ! Write the closing tags
542              WRITE(PVD_UNIT,"(3X,A)")'</Collection>'
543              WRITE(PVD_UNIT,"(A)")'</VTKFile>'
544     
545              CLOSE(PVD_UNIT)
546           ENDIF
547     ! Identify that the files has been created and opened for next pass
548           FIRST_PASS = .FALSE.
549     
550           CALL FINL_ERR_MSG
551     
552     ! Return to the calling routine
553           RETURN
554     
555           END SUBROUTINE ADD_VTP_TO_PVD
556     
557     
558     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
559     !                                                                      C
560     !  Module name: WRITE_VTP_FILE                                         C
561     !  Purpose: Writes particles data in VTK format (Polydata VTP)         C
562     !           Binary format                                              C
563     !                                                                      C
564     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
565     !  Reviewer:                                          Date:            C
566     !                                                                      C
567     !  Revision Number #                                  Date: ##-###-##  C
568     !  Author: #                                                           C
569     !  Purpose: #                                                          C
570     !                                                                      C
571     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
572           SUBROUTINE WRITE_VTP_FILE(LCV,MODE)
573     
574           USE vtk, only: DIMENSION_VTK, VTK_DEFINED, FRAME
575           USE vtk, only: VTK_REGION,VTK_DEFINED,VTK_DATA
576           USE vtk, only: VTK_PART_DIAMETER,VTK_PART_VEL,VTK_PART_USR_VAR,VTK_PART_TEMP
577           USE vtk, only: VTK_PART_ANGULAR_VEL,VTK_PART_ORIENTATION
578           USE vtk, only: VTK_PART_X_S, VTK_PART_COHESION
579           USE vtk, only: TIME_DEPENDENT_FILENAME,VTU_FRAME_UNIT,VTU_FRAME_FILENAME
580           USE vtk, only: VTK_DBG_FILE
581           USE output, only: FULL_LOG
582           use des_thermo, only: DES_T_s
583           IMPLICIT NONE
584           INTEGER :: L,N,LCV
585     
586           INTEGER :: PASS
587           INTEGER :: WRITE_HEADER = 1
588           INTEGER :: WRITE_DATA   = 2
589           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
590                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
591     
592           VTK_REGION = LCV
593     ! There is nothing to write if we are not in a defined vtk region
594           IF(.NOT.VTK_DEFINED(VTK_REGION)) RETURN
595     
596           IF(VTK_DATA(LCV)/='P') RETURN
597           IF(MODE==0.AND.(VTK_DBG_FILE(LCV))) RETURN
598           IF(MODE==1.AND.(.NOT.VTK_DBG_FILE(LCV))) RETURN
599     
600           CALL SETUP_VTK_REGION_PARTICLES
601     
602           CALL OPEN_VTP_FILE_BIN(MODE)
603     
604     ! Only open pvd file when there are particles in vtk region
605           IF(GLOBAL_CNT>0.AND.MODE==0) CALL OPEN_PVD_FILE
606     
607     ! First pass write the data header.
608     ! Second pass writes the data (appended binary format).
609     
610           DO PASS=WRITE_HEADER,WRITE_DATA
611     
612     
613              CALL WRITE_GEOMETRY_IN_VTP_BIN(PASS)
614     
615              IF(VTK_PART_DIAMETER(VTK_REGION)) &
616                 CALL WRITE_SCALAR_IN_VTP_BIN('Diameter',2.0D0*DES_RADIUS,PASS)
617     
618              IF(VTK_PART_VEL(VTK_REGION)) &
619                 CALL WRITE_VECTOR_IN_VTP_BIN('Velocity',DES_VEL_NEW,PASS)
620     
621              IF(VTK_PART_ANGULAR_VEL(VTK_REGION)) &
622                 CALL WRITE_VECTOR_IN_VTP_BIN('Angular_velocity', OMEGA_NEW,PASS)
623     
624              IF(PARTICLE_ORIENTATION) THEN
625                 IF(VTK_PART_ORIENTATION(VTK_REGION)) &
626                    CALL WRITE_VECTOR_IN_VTP_BIN('Orientation', ORIENTATION,PASS)
627              ENDIF
628     
629              DO N=1, DES_USR_VAR_SIZE
630                 IF(VTK_PART_USR_VAR(VTK_REGION,N)) &
631                   CALL WRITE_SCALAR_IN_VTP_BIN('User Defined Var '//trim(iVal(N)),DES_USR_VAR(N,:),PASS)
632              ENDDO
633     
634              IF(ENERGY_EQ.AND.VTK_PART_TEMP(VTK_REGION)) &
635                CALL WRITE_SCALAR_IN_VTP_BIN('Temperature', DES_T_s,PASS)
636     
637              IF(ANY_SPECIES_EQ) THEN
638                 DO N=1, DIMENSION_N_S
639                    IF(VTK_PART_X_s(VTK_REGION,N)) &
640                      CALL WRITE_SCALAR_IN_VTP_BIN(trim(iVar('X_s',N)), DES_X_s(:,N),PASS)
641                 ENDDO
642              ENDIF
643     
644           IF(USE_COHESION.AND.VTK_PART_COHESION(VTK_REGION)) &
645              CALL WRITE_SCALAR_IN_VTP_BIN('CohesiveForce', PostCohesive,PASS)
646     
647           ENDDO ! PASS LOOP, EITHER HEADER OR DATA
648     
649     
650           CALL CLOSE_VTP_FILE_BIN(MODE)
651     
652     ! Only update pvd file when there are particles in vtk region
653           IF(GLOBAL_CNT>0.AND.MODE==0) CALL UPDATE_AND_CLOSE_PVD_FILE
654     
655     #ifdef MPI
656           call MPI_barrier(MPI_COMM_WORLD,mpierr)
657     #endif
658     
659     ! Update Frames
660           IF (myPE == PE_IO.AND.TIME_DEPENDENT_FILENAME) THEN
661              OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_FRAME_UNIT, FILE = TRIM(VTU_FRAME_FILENAME))
662              DO L = 1, DIMENSION_VTK
663                 IF(VTK_DEFINED(L)) WRITE(VTU_FRAME_UNIT,*) L,FRAME(L)
664              ENDDO
665              CLOSE(VTU_FRAME_UNIT)
666           ENDIF
667     
668          IF (FULL_LOG.AND.myPE == PE_IO) WRITE(*,20)' DONE.'
669     
670     20    FORMAT(A,1X/)
671           RETURN
672     
673           END SUBROUTINE WRITE_VTP_FILE
674     
675     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
676     !                                                                      C
677     !  Module name: OPEN_VTP_FILE                                          C
678     !  Purpose: Open a vtp file and writes the header                      C
679     !           Binary format                                              C
680     !                                                                      C
681     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
682     !  Reviewer:                                          Date:            C
683     !                                                                      C
684     !  Revision Number #                                  Date: ##-###-##  C
685     !  Author: #                                                           C
686     !  Purpose: #                                                          C
687     !                                                                      C
688     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
689       SUBROUTINE OPEN_VTP_FILE_BIN(MODE)
690     
691           USE run, only: TIME
692           USE output, only: FULL_LOG
693           USE vtk, only: TIME_DEPENDENT_FILENAME, VTU_FRAME_FILENAME, VTU_FRAME_UNIT
694           USE vtk, only: RESET_FRAME_AT_TIME_ZERO,PVTU_FILENAME,PVTU_UNIT,BUFFER,END_REC
695           USE vtk, only: DIMENSION_VTK, VTK_DEFINED,  FRAME,  VTK_REGION
696           USE vtk, only: VTU_FILENAME, VTK_FILEBASE, VTU_DIR, VTU_UNIT
697           USE param1, only: ZERO
698     
699           IMPLICIT NONE
700           LOGICAL :: VTU_FRAME_FILE_EXISTS, NEED_TO_WRITE_VTP
701           INTEGER :: ISTAT,BUFF1,BUFF2,L
702           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
703                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
704     
705     
706           IF(BDIST_IO) THEN
707              NEED_TO_WRITE_VTP = (LOCAL_CNT>0)
708           ELSE
709              NEED_TO_WRITE_VTP = (MyPE==0.AND.GLOBAL_CNT>0)
710           ENDIF
711     
712     ! Only open the file from head node when not using distributed I/O
713           IF (myPE /= PE_IO.AND.(.NOT.BDIST_IO)) RETURN
714     
715           IF(TIME_DEPENDENT_FILENAME) THEN
716              INQUIRE(FILE=VTU_FRAME_FILENAME,EXIST=VTU_FRAME_FILE_EXISTS)
717              IF(VTU_FRAME_FILE_EXISTS) THEN
718                 OPEN(CONVERT='BIG_ENDIAN',UNIT = VTU_FRAME_UNIT, FILE = TRIM(VTU_FRAME_FILENAME))
719                 DO L = 1, DIMENSION_VTK
720                    IF(VTK_DEFINED(L)) THEN
721                       READ(VTU_FRAME_UNIT,*)BUFF1,BUFF2
722                       FRAME(L)=BUFF2
723                    ENDIF
724                 ENDDO
725                 CLOSE(VTU_FRAME_UNIT)
726              ENDIF
727              IF(RESET_FRAME_AT_TIME_ZERO.AND.TIME==ZERO) THEN
728                 DO L = 1, DIMENSION_VTK
729                    IF(L==VTK_REGION) FRAME(L)=-1
730                 ENDDO
731              ENDIF
732              DO L = 1, DIMENSION_VTK
733                 IF(L==VTK_REGION) FRAME(L) = FRAME(L) + 1
734              ENDDO
735           ENDIF
736     
737     ! For distributed I/O, define the file name for each processor that owns particles
738           IF (BDIST_IO) THEN
739              IF (LOCAL_CNT>0) THEN
740                 IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
741                    WRITE(VTU_FILENAME,20) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),MYPE
742                 ELSE
743                    WRITE(VTU_FILENAME,25) TRIM(VTK_FILEBASE(VTK_REGION)),MYPE
744                 ENDIF
745              ENDIF
746           ELSE
747              IF(MYPE.EQ.PE_IO) THEN
748                 IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
749                    WRITE(VTU_FILENAME,30) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION)
750                 ELSE
751                    WRITE(VTU_FILENAME,35) TRIM(VTK_FILEBASE(VTK_REGION))
752                 ENDIF
753              END IF
754           END IF
755     
756     ! Add the VTU directory path if necessary
757     
758           IF (NEED_TO_WRITE_VTP) THEN
759              IF(TRIM(VTU_DIR)/='.') VTU_FILENAME='./'//TRIM(VTU_DIR)//'/'//VTU_FILENAME
760           ENDIF
761     
762     ! Echo
763           IF (FULL_LOG) THEN
764              IF (.NOT.BDIST_IO) THEN
765                 WRITE(*,10,ADVANCE='NO')' WRITING VTP FILE : ', TRIM(VTU_FILENAME),' .'
766              ELSE
767                 IF(myPE==PE_IO) WRITE(*,15,ADVANCE='NO')' EACH PROCESOR IS WRITING ITS OWN VTP FILE.'
768              ENDIF
769           ENDIF
770     
771     ! Open File
772     
773           IF (NEED_TO_WRITE_VTP) THEN
774     
775              VTU_UNIT = 678
776              OPEN(CONVERT='BIG_ENDIAN',UNIT     = VTU_UNIT,           &
777                   FILE     = TRIM(VTU_FILENAME), &
778                   FORM     = 'UNFORMATTED',      &  ! works with gfortran 4.3.4 and ifort 10.1 but may not be supported by all compilers
779                                                     ! use 'BINARY' if 'UNFORMATTED' is not supported
780                   ACCESS   = 'STREAM',           &  ! works with gfortran 4.3.4 and ifort 10.1 but may not be supported by all compilers
781                                                     ! use 'SEQUENTIAL' if 'STREAM' is not supported
782                   ACTION   = 'WRITE',            &
783                   IOSTAT=ISTAT)
784     
785     
786              IF (ISTAT /= 0) THEN
787                 IF(DMP_LOG) WRITE(UNIT_LOG, 1001) VTU_FILENAME, VTU_UNIT,VTU_DIR
788                 IF(FULL_LOG.AND.myPE == PE_IO) WRITE(*, 1001) VTU_FILENAME, VTU_UNIT,VTU_DIR
789                 CALL MFIX_EXIT(myPE)
790              ENDIF
791     
792     
793     1001 FORMAT(/1X,70('*')//, ' From: OPEN_VTP_FILE',/,' Message: ',          &
794                 'Error opening vtp file. Terminating run.',/10X,          &
795                 'File name:  ',A,/10X,                                         &
796                 'VTP_UNIT :  ',i4, /10X,                                       &
797                 'PLEASE VERIFY THAT VTU_DIR EXISTS: ', A, &
798                 /1X,70('*')/)
799     
800     
801     ! Write file Header
802              BUFFER='<?xml version="1.0"?>'
803              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
804     
805              WRITE(BUFFER,*)'<!-- Time =',TIME,' sec. -->'
806              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
807     
808              BUFFER='<VTKFile type="PolyData" version="0.1" byte_order="BigEndian">'
809              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
810     
811              BUFFER='  <PolyData>'
812              WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
813     
814     
815           ENDIF
816     ! For distributed I/O, open .p))vtp file that combines all *.vtp files for a given FRAME
817     ! this is a simple ASCII file
818     
819           IF (myPE == PE_IO.AND.BDIST_IO.AND.GLOBAL_CNT>0) THEN
820     
821              IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
822                 WRITE(PVTU_FILENAME,40) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION)
823              ELSE
824                 WRITE(PVTU_FILENAME,45) TRIM(VTK_FILEBASE(VTK_REGION))
825              ENDIF
826     
827              IF(TRIM(VTU_DIR)/='.') PVTU_FILENAME='./'//TRIM(VTU_DIR)//'/'//PVTU_FILENAME
828     
829              OPEN(CONVERT='BIG_ENDIAN',UNIT = PVTU_UNIT, FILE = TRIM(PVTU_FILENAME))
830     
831              WRITE(PVTU_UNIT,100) '<?xml version="1.0"?>'
832              WRITE(PVTU_UNIT,110) '<!-- Time =',TIME,' sec. -->'
833              WRITE(PVTU_UNIT,120) '<VTKFile type="PPolyData"',&
834                       ' version="0.1" byte_order="BigEndian">'
835     
836              WRITE(PVTU_UNIT,100) '  <PPolyData GhostLevel="0">'
837              WRITE(PVTU_UNIT,100) '      <PPoints>'
838              WRITE(PVTU_UNIT,100) '        <PDataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
839                   &format="appended" offset=" 0" />'
840              WRITE(PVTU_UNIT,100) '      </PPoints>'
841              WRITE(PVTU_UNIT,100) ''
842              WRITE(PVTU_UNIT,100) '      <PPointData Scalars="Diameter" Vectors="Velocity">'
843     
844           ENDIF
845     
846     100   FORMAT(A)
847     110   FORMAT(A,E14.7,A)
848     120   FORMAT(A,A)
849     10    FORMAT(/1X,3A)
850     15    FORMAT(/1X,A)
851     20    FORMAT(A,"_",I4.4,"_",I5.5,".vtp")
852     25    FORMAT(A,"_",I5.5,".vtp")
853     30    FORMAT(A,"_",I4.4,".vtp")
854     35    FORMAT(A,".vtp")
855     40    FORMAT(A,"_",I4.4,".pvtp")
856     45    FORMAT(A,".pvtp")
857     
858           RETURN
859     
860           END SUBROUTINE OPEN_VTP_FILE_BIN
861     
862     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
863     !                                                                      C
864     !  Module name: WRITE_GEOMETRY_IN_VTP_BIN                              C
865     !  Purpose: Write Geometry and connectivity in a vtu file              C
866     !           Binary format                                              C
867     !                                                                      C
868     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
869     !  Reviewer:                                          Date:            C
870     !                                                                      C
871     !  Revision Number #                                  Date: ##-###-##  C
872     !  Author: #                                                           C
873     !  Purpose: #                                                          C
874     !                                                                      C
875     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
876       SUBROUTINE WRITE_GEOMETRY_IN_VTP_BIN(PASS)
877     
878           USE vtk, only: NUMBER_OF_POINTS,BUFFER, VTU_UNIT,END_REC,VTU_OFFSET,BELONGS_TO_VTK_SUBDOMAIN
879     
880           IMPLICIT NONE
881     
882           REAL(c_float) :: float
883           INTEGER(c_int) :: int
884     
885           INTEGER ::     nbytes_vector
886           INTEGER ::     offset_xyz
887     
888           INTEGER :: PASS
889           INTEGER :: WRITE_HEADER = 1
890           INTEGER :: WRITE_DATA   = 2
891     
892           DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)  ! local
893           DOUBLE PRECISION, ALLOCATABLE :: gtemp_array(:,:)  ! global
894     
895           INTEGER :: LB, UB
896           INTEGER :: PC, LC1, LC2
897     
898     ! Loop through all particles and kee a list of particles belonging to a VTK region
899     
900     ! Since the data is appended (i.e., written after all tags), the
901     ! offset, in number of bytes must be specified.  The offset includes
902     ! the size of the data for each field, plus the size of the integer
903     ! that stores the number of bytes.  this is why the offset of a field
904     ! equals the offset of the previous field plus sizeof(int) plus the
905     ! number of bytes of the field.
906     
907     ! Next, the actual data is written for the geometry (PASS=WRITE_DATA)
908     ! The DATA is converted to single precision to save memory.
909     
910           IF (.NOT.BDIST_IO) THEN
911     ! The number of points in the pvd file is the global number of particles
912     ! computed from SETUP_VTK_REGION_PARTICLES
913     
914              NUMBER_OF_POINTS = GLOBAL_CNT
915     
916     ! Number of bytes of position field (vector,3 components)
917              nbytes_vector       = NUMBER_OF_POINTS * 3 * sizeof(float)
918     
919     ! Offset of each field
920              offset_xyz = 0
921     
922              IF(PASS==WRITE_HEADER) THEN
923                 IF(myPE == PE_IO) THEN
924     
925                    WRITE(BUFFER,*)'    <Piece NumberOfPoints="',NUMBER_OF_POINTS, &
926                          '"  NumberOfVerts="0" NumberOfLines ="0" NumberOfStrips="0" NumberOfPolys="0" >'
927                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
928     
929                    WRITE(BUFFER,*)'      <Points>'
930                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
931     
932                    WRITE(BUFFER,*)'        <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
933                                            &format="appended" offset="',offset_xyz,'" />'
934                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
935     
936                    WRITE(BUFFER,*)'      </Points>'
937                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
938     
939                    WRITE(BUFFER,*)'<PointData Scalars="Diameter" Vectors="Velocity"> '!preparing pointData
940                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
941     
942     ! calculate offset for next field
943                    VTU_offset = offset_xyz + sizeof(int) + nbytes_vector
944     
945                 ENDIF
946     
947              ELSEIF(PASS==WRITE_DATA) THEN
948     
949                 IF(myPE == PE_IO) THEN
950     
951                    WRITE(BUFFER,*)'      </PointData>'
952                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
953     
954                    WRITE(BUFFER,*)'      <Verts> </Verts>'
955                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
956     
957                    WRITE(BUFFER,*)'      <Lines> </Lines>'
958                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
959     
960                    WRITE(BUFFER,*)'      <Strips> </Strips>'
961                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
962     
963                    WRITE(BUFFER,*)'      <Polys> </Polys>'
964                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
965     
966                    WRITE(BUFFER,*)'    </Piece>'
967                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
968     
969                    WRITE(BUFFER,*)'  </PolyData>'
970                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
971     
972                    WRITE(BUFFER,*)'  <AppendedData encoding="raw">'
973                    WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
974     
975     
976     ! Starting raw binary data with an underscore
977     
978                    WRITE(BUFFER,*)'_'
979                    WRITE(VTU_UNIT)TRIM(BUFFER)
980     
981     ! Number of bytes for X,Y,Z coordinates
982                 WRITE(VTU_UNIT) nbytes_vector
983     
984     
985              ENDIF
986     
987              LB = LBOUND(DES_POS_NEW,2) ! This should always be 1
988              UB = UBOUND(DES_POS_NEW,2) ! This should always be 2
989     
990              ALLOCATE (dProcBuf(LOCAL_CNT) )
991              ALLOCATE (dRootBuf(GLOBAL_CNT))
992              ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
993              ALLOCATE (gtemp_array((UB-LB)+1,GLOBAL_CNT))
994     
995     ! Pack particle coordinates in a temporary local array
996              PC = 0
997              DO LC1 = 1, MAX_PIP
998                 IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
999                    PC =PC + 1
1000                    DO LC2=LB, UB
1001                       ltemp_array(LC2,PC) = DES_POS_NEW(LC1,LC2)
1002                    ENDDO
1003                 ENDIF
1004                 IF(PC==LOCAL_CNT) EXIT
1005              ENDDO
1006     
1007     ! For each coordinate (x,y, and z), gather the local list to global temporary array
1008              DO LC1 = LB, UB
1009                 dprocbuf(1:LOCAL_CNT)=ltemp_array(LC1,1:LOCAL_CNT)
1010                 CALL desmpi_gatherv(ptype=2)
1011                 gtemp_array(LC1,:) = drootbuf(:)
1012              ENDDO
1013     
1014     ! Write the list of coordinates
1015              IF(myPE == PE_IO) THEN
1016                 DO LC1=1, GLOBAL_CNT
1017                    DO LC2=LB, UB
1018                       WRITE(VTU_UNIT)  real(gtemp_array(LC2,LC1))
1019                    ENDDO
1020                 ENDDO
1021              ENDIF
1022     
1023              deallocate (dProcBuf, dRootBuf, ltemp_array,gtemp_array)
1024     
1025     
1026              ENDIF
1027     
1028     
1029           ELSEIF(BDIST_IO.AND.LOCAL_CNT>0) THEN
1030     
1031              IF(LOCAL_CNT==0) RETURN
1032     ! The number of points in the pvd file is the local number of particles
1033     ! computed from SETUP_VTK_REGION_PARTICLES
1034     
1035              NUMBER_OF_POINTS = LOCAL_CNT
1036     
1037     ! Number of bytes of position field (vector,3 components)
1038              nbytes_vector       = NUMBER_OF_POINTS * 3 * sizeof(float)
1039     
1040     ! Offset of each field
1041              offset_xyz = 0
1042     
1043              IF(PASS==WRITE_HEADER) THEN
1044     
1045                 WRITE(BUFFER,*)'    <Piece NumberOfPoints="',NUMBER_OF_POINTS, &
1046                       '"  NumberOfVerts="0" NumberOfLines ="0" NumberOfStrips="0" NumberOfPolys="0" >'
1047                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1048     
1049                 WRITE(BUFFER,*)'      <Points>'
1050                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1051     
1052                 WRITE(BUFFER,*)'        <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" &
1053                                         &format="appended" offset="',offset_xyz,'" />'
1054                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1055     
1056                 WRITE(BUFFER,*)'      </Points>'
1057                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1058     
1059                 WRITE(BUFFER,*)'<PointData Scalars="Diameter" Vectors="Velocity"> '!preparing pointData
1060                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1061     
1062     ! calculate offset for next field
1063                 VTU_offset = offset_xyz + sizeof(int) + nbytes_vector
1064     
1065     
1066              ELSEIF(PASS==WRITE_DATA) THEN
1067     
1068                 WRITE(BUFFER,*)'      </PointData>'
1069                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1070     
1071                 WRITE(BUFFER,*)'      <Verts> </Verts>'
1072                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1073     
1074                 WRITE(BUFFER,*)'      <Lines> </Lines>'
1075                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1076     
1077                 WRITE(BUFFER,*)'      <Strips> </Strips>'
1078                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1079     
1080                 WRITE(BUFFER,*)'      <Polys> </Polys>'
1081                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1082     
1083                 WRITE(BUFFER,*)'    </Piece>'
1084                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1085     
1086                 WRITE(BUFFER,*)'  </PolyData>'
1087                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1088     
1089                 WRITE(BUFFER,*)'  <AppendedData encoding="raw">'
1090                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1091     
1092     ! Starting raw binary data with an underscore
1093     
1094                 WRITE(BUFFER,*)'_'
1095                 WRITE(VTU_UNIT)TRIM(BUFFER)
1096     
1097     ! Number of bytes for X,Y,Z coordinates
1098                 WRITE(VTU_UNIT) nbytes_vector
1099     
1100                 LB = LBOUND(DES_POS_NEW,2) ! This should always be 1
1101                 UB = UBOUND(DES_POS_NEW,2) ! This should always be 2
1102     
1103                 ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
1104     
1105     ! Pack particle coordinates in a temporary local array
1106                 PC = 0
1107                 DO LC1 = 1, MAX_PIP
1108                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1109                       PC =PC + 1
1110                       DO LC2=LB, UB
1111                          ltemp_array(LC2,PC) = DES_POS_NEW(LC1,LC2)
1112                       ENDDO
1113                    ENDIF
1114                    IF(PC==LOCAL_CNT) EXIT
1115                 ENDDO
1116     
1117     
1118     ! Write the list of coordinates
1119                 DO LC1=1, LOCAL_CNT
1120                    DO LC2=LB, UB
1121                       WRITE(VTU_UNIT)  real(ltemp_array(LC2,LC1))
1122                    ENDDO
1123                 ENDDO
1124     
1125                 deallocate (ltemp_array)
1126     
1127              ENDIF
1128     
1129           ENDIF
1130     
1131           RETURN
1132     
1133           END SUBROUTINE WRITE_GEOMETRY_IN_VTP_BIN
1134     
1135     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1136     !                                                                      C
1137     !  Module name: WRITE_SCALAR_IN_VTP_BIN                                C
1138     !  Purpose: Write Scalar variable in a vtp file                        C
1139     !           Binary format                                              C
1140     !                                                                      C
1141     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1142     !  Reviewer:                                          Date:            C
1143     !                                                                      C
1144     !  Revision Number #                                  Date: ##-###-##  C
1145     !  Author: #                                                           C
1146     !  Purpose: #                                                          C
1147     !                                                                      C
1148     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1149       SUBROUTINE WRITE_SCALAR_IN_VTP_BIN(VAR_NAME,VAR,PASS)
1150     
1151           USE vtk, only: BUFFER,VTU_OFFSET,VTU_UNIT,PVTU_UNIT
1152           USE vtk, only: END_REC,BELONGS_TO_VTK_SUBDOMAIN
1153           USE output, only: FULL_LOG
1154     
1155           IMPLICIT NONE
1156           INTEGER :: I,LC1,PC
1157     
1158           CHARACTER (*) :: VAR_NAME
1159           DOUBLE PRECISION, INTENT(in) :: VAR(:)
1160     
1161           REAL(c_float) :: float
1162     
1163           INTEGER :: nbytes_scalar
1164     
1165           INTEGER :: PASS
1166           INTEGER :: WRITE_HEADER = 1
1167           INTEGER :: WRITE_DATA   = 2
1168     
1169           IF (.NOT.BDIST_IO) THEN
1170     
1171     ! Number of bytes for each scalar field
1172              nbytes_scalar = GLOBAL_CNT * sizeof(float)
1173     
1174              IF(PASS==WRITE_HEADER) THEN
1175     
1176     ! Remove possible white space with underscore
1177                 DO I = 1,LEN_TRIM(VAR_NAME)
1178                    IF(VAR_NAME(I:I) == ' ') VAR_NAME(I:I) = '_'
1179                 ENDDO
1180     
1181     ! For each scalar, write a tag, with corresponding offset
1182                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1183                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1184                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1185     
1186     ! Prepare the offset for the next field
1187                 VTU_offset = VTU_offset + sizeof(float) + nbytes_scalar
1188     
1189     
1190              ELSEIF(PASS==WRITE_DATA) THEN
1191     
1192                allocate (dProcBuf(LOCAL_CNT) )
1193                allocate (dRootBuf(GLOBAL_CNT))
1194     
1195     ! Pack scalar list in a local buffer before gathering to root
1196                 PC = 0
1197                 DO LC1 = 1, MAX_PIP
1198                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1199                       PC =PC + 1
1200                       dProcBuf(PC) = VAR(LC1)
1201                    ENDIF
1202                    IF(PC==LOCAL_CNT) EXIT
1203                 ENDDO
1204     
1205     ! Gather local buffer to root
1206              CALL desmpi_gatherv(ptype=2)
1207     
1208     ! Write the data, always preceded by its size in number of bytes
1209     ! Write root buffer to file
1210              WRITE(VTU_UNIT) nbytes_scalar
1211     
1212              IF(myPE == PE_IO) THEN
1213                 DO LC1=1, GLOBAL_CNT
1214                    WRITE(VTU_UNIT)  real(drootBuf(LC1))
1215                 ENDDO
1216              ENDIF
1217     
1218              deallocate (dProcBuf, dRootBuf)
1219     
1220     
1221              ENDIF
1222     
1223     
1224           ELSEIF(BDIST_IO.AND.LOCAL_CNT>0) THEN
1225     
1226     ! Number of bytes for each scalar field
1227              nbytes_scalar = LOCAL_CNT * sizeof(float)
1228     
1229     ! Remove possible white space with underscore
1230              DO I = 1,LEN_TRIM(VAR_NAME)
1231                 IF(VAR_NAME(I:I) == ' ') VAR_NAME(I:I) = '_'
1232              ENDDO
1233     
1234              IF(PASS==WRITE_HEADER) THEN
1235     
1236     ! For each scalar, write a tag, with corresponding offset
1237                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1238                      TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1239                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1240     
1241     ! Prepare the offset for the next field
1242                 VTU_offset = VTU_offset + sizeof(float) + nbytes_scalar
1243     
1244     
1245              ELSEIF(PASS==WRITE_DATA) THEN
1246     
1247                 allocate (dProcBuf(LOCAL_CNT) )
1248     
1249     ! Pack scalar list in a local buffer before writing in file
1250                 PC = 0
1251                 DO LC1 = 1, MAX_PIP
1252                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1253                       PC =PC + 1
1254                       dProcBuf(PC) = VAR(LC1)
1255                    ENDIF
1256                    IF(PC==LOCAL_CNT) EXIT
1257                 ENDDO
1258     
1259     ! Write the data, always preceded by its size in number of bytes
1260     ! Write root buffer to file
1261                 WRITE(VTU_UNIT) nbytes_scalar
1262     
1263                 DO LC1=1, LOCAL_CNT
1264                    WRITE(VTU_UNIT)  real(dProcBuf(LC1))
1265                 ENDDO
1266     
1267                 deallocate (dProcBuf)
1268     
1269                 IF (myPE == PE_IO) THEN       ! Update pvtu file with variable name
1270                    WRITE(PVTU_UNIT,90) '        <PointArray type="Float32" Name="', &
1271                         TRIM(VAR_NAME),'" format="appended" offset="',VTU_offset,'" />'
1272                 ENDIF
1273     
1274     
1275              ENDIF
1276     
1277     
1278           ENDIF
1279     
1280     
1281           IF (PASS==WRITE_DATA.AND.FULL_LOG.AND.myPE == PE_IO) WRITE(*,10,ADVANCE='NO')'.'
1282     
1283     10    FORMAT(A)
1284     90    FORMAT(A,A,A,I12,A)
1285     
1286           RETURN
1287     
1288           END SUBROUTINE WRITE_SCALAR_IN_VTP_BIN
1289     
1290     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1291     !                                                                      C
1292     !  Module name: WRITE_VECTOR_IN_VTP                                    C
1293     !  Purpose: Write Vector variable in a vtp file                        C
1294     !           Binary format                                              C
1295     !                                                                      C
1296     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1297     !  Reviewer:                                          Date:            C
1298     !                                                                      C
1299     !  Revision Number #                                  Date: ##-###-##  C
1300     !  Author: #                                                           C
1301     !  Purpose: #                                                          C
1302     !                                                                      C
1303     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1304       SUBROUTINE WRITE_VECTOR_IN_VTP_BIN(VAR_NAME,VAR,PASS)
1305     
1306           USE vtk, only: BUFFER,VTU_OFFSET,VTU_UNIT,PVTU_UNIT
1307           USE vtk, only: END_REC,BELONGS_TO_VTK_SUBDOMAIN
1308           USE output, only: FULL_LOG
1309     
1310           IMPLICIT NONE
1311     
1312           CHARACTER (*) :: VAR_NAME
1313           DOUBLE PRECISION, INTENT(in) :: VAR(:,:)
1314     
1315           REAL(c_float) :: float
1316     
1317           INTEGER :: nbytes_vector
1318     
1319           INTEGER :: PASS
1320           INTEGER :: WRITE_HEADER = 1
1321           INTEGER :: WRITE_DATA   = 2
1322     
1323           DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)  ! local
1324           DOUBLE PRECISION, ALLOCATABLE :: gtemp_array(:,:)  ! global
1325     
1326           INTEGER :: LB, UB
1327           INTEGER :: PC, LC1, LC2
1328     
1329           IF (.NOT.BDIST_IO) THEN
1330     
1331     ! Number of bytes for each vector field
1332              nbytes_vector = GLOBAL_CNT * 3 * sizeof(float)
1333     
1334              IF(PASS==WRITE_HEADER) THEN
1335     ! For each vector, write a tag, with corresponding offset
1336     
1337                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1338                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1339                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1340     
1341     ! Prepare the offset for the next field
1342                 VTU_offset = VTU_offset + sizeof(float) + nbytes_vector
1343     
1344     
1345              ELSEIF(PASS==WRITE_DATA) THEN
1346     
1347                 LB = LBOUND(VAR,2) ! This should always be 1
1348                 UB = UBOUND(VAR,2) ! This should always be 2
1349     
1350                 ALLOCATE (dProcBuf(LOCAL_CNT) )
1351                 ALLOCATE (dRootBuf(GLOBAL_CNT))
1352                 ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
1353                 ALLOCATE (gtemp_array((UB-LB)+1,GLOBAL_CNT))
1354     
1355     ! For each vector component, pack component list in a local array
1356                 PC = 0
1357                 DO LC1 = 1, MAX_PIP
1358                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1359                       PC =PC + 1
1360                       DO LC2=LB, UB
1361                          ltemp_array(LC2,PC) = VAR(LC1,LC2)
1362                       ENDDO
1363                    ENDIF
1364                    IF(PC==LOCAL_CNT) EXIT
1365                 ENDDO
1366     
1367     
1368     ! For each component, gather the local list to global temporary array
1369              DO LC1 = LB, UB
1370                 dprocbuf(1:LOCAL_CNT)=ltemp_array(LC1,1:LOCAL_CNT)
1371                 CALL desmpi_gatherv(ptype=2)
1372                 gtemp_array(LC1,:) = drootbuf(:)
1373              ENDDO
1374     
1375     ! Write the data, always preceded by its size in number of bytes
1376              IF(myPE == PE_IO) THEN
1377                 WRITE(VTU_UNIT) nbytes_vector
1378                 DO LC1=1, GLOBAL_CNT
1379                    DO LC2=LB, UB
1380                       WRITE(VTU_UNIT)  real(gtemp_array(LC2,LC1))
1381                    ENDDO
1382                 ENDDO
1383              ENDIF
1384     
1385              deallocate (dProcBuf, dRootBuf, ltemp_array,gtemp_array)
1386     
1387     
1388     
1389              ENDIF
1390     
1391     
1392           ELSEIF(BDIST_IO.AND.LOCAL_CNT>0) THEN
1393     
1394     ! Number of bytes for each vector field
1395              nbytes_vector = LOCAL_CNT * 3 * sizeof(float)
1396     
1397              IF(PASS==WRITE_HEADER) THEN
1398     ! For each vector, write a tag, with corresponding offset
1399     
1400                 WRITE(BUFFER,90)'        <DataArray type="Float32" Name="', &
1401                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1402                 WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1403     
1404     ! Prepare the offset for the next field
1405                 VTU_offset = VTU_offset + sizeof(float) + nbytes_vector
1406     
1407     
1408              ELSEIF(PASS==WRITE_DATA) THEN
1409     
1410                 LB = LBOUND(VAR,1) ! This should always be 1
1411                 UB = UBOUND(VAR,1) ! This should always be 2
1412     
1413                 ALLOCATE (ltemp_array((UB-LB)+1,LOCAL_CNT))
1414     
1415     ! For each vector component, pack component list in a local array
1416                 PC = 0
1417                 DO LC1 = 1, MAX_PIP
1418                    IF(BELONGS_TO_VTK_SUBDOMAIN(LC1)) THEN
1419                       PC =PC + 1
1420                       DO LC2=LB, UB
1421                          ltemp_array(LC2,PC) = VAR(LC2,LC1)
1422                       ENDDO
1423                    ENDIF
1424                    IF(PC==LOCAL_CNT) EXIT
1425                 ENDDO
1426     
1427     
1428     ! Write the data, always preceded by its size in number of bytes
1429                 WRITE(VTU_UNIT) nbytes_vector
1430                 DO LC1=1, LOCAL_CNT
1431                    DO LC2=LB, UB
1432                       WRITE(VTU_UNIT)  real(ltemp_array(LC2,LC1))
1433                    ENDDO
1434                 ENDDO
1435     
1436                 deallocate (ltemp_array)
1437     
1438     
1439                 IF (myPE == PE_IO) THEN       ! Update pvtu file with variable name
1440                    WRITE(PVTU_UNIT,90)'        <PointArray type="Float32" Name="', &
1441                      TRIM(VAR_NAME),'"  NumberOfComponents="3" format="appended" offset="',VTU_offset,'" />'
1442                 ENDIF
1443     
1444     
1445              ENDIF
1446     
1447           ENDIF
1448     
1449     
1450           IF (PASS==WRITE_DATA.AND.FULL_LOG.AND.myPE == PE_IO) WRITE(*,10,ADVANCE='NO')'.'
1451     
1452     10    FORMAT(A)
1453     90    FORMAT(A,A,A,I12,A)
1454     
1455           RETURN
1456     
1457           END SUBROUTINE WRITE_VECTOR_IN_VTP_BIN
1458     
1459     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1460     !                                                                      C
1461     !  Module name: CLOSE_VTP_FILE_BIN                                     C
1462     !  Purpose: Close a vtp file                                           C
1463     !           Binary format                                              C
1464     !                                                                      C
1465     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1466     !  Reviewer:                                          Date:            C
1467     !                                                                      C
1468     !  Revision Number #                                  Date: ##-###-##  C
1469     !  Author: #                                                           C
1470     !  Purpose: #                                                          C
1471     !                                                                      C
1472     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1473       SUBROUTINE CLOSE_VTP_FILE_BIN(MODE)
1474     
1475           USE vtk, only: BUFFER,VTU_UNIT,END_REC,PVTU_UNIT,TIME_DEPENDENT_FILENAME
1476           USE vtk, only: VTK_REGION,VTK_FILEBASE,FRAME
1477     
1478           IMPLICIT NONE
1479     
1480           INTEGER:: N
1481           CHARACTER (LEN=32)  :: VTU_NAME
1482           INTEGER, DIMENSION(0:numPEs-1) :: ALL_PART_CNT
1483           INTEGER :: IERR
1484           INTEGER :: MODE   ! MODE = 0 : Write regular VTK region file
1485                             ! MODE = 1 : Write debug   VTK region file (VTK_DBG_FILE = .TRUE.)
1486     
1487     
1488           IF((myPE == PE_IO.AND.(.NOT.BDIST_IO)).OR.(BDIST_IO.AND.LOCAL_CNT>0)) THEN
1489     
1490     ! Write last tags and close the vtp file
1491           WRITE(BUFFER,110)'  </AppendedData>'
1492           WRITE(VTU_UNIT)END_REC//TRIM(BUFFER)//END_REC
1493     
1494           WRITE(BUFFER,110)'</VTKFile>'
1495           WRITE(VTU_UNIT)TRIM(BUFFER)//END_REC
1496     
1497           CLOSE(VTU_UNIT)
1498     
1499           ENDIF
1500     
1501     ! Update pvtu file and close
1502     
1503           IF(BDIST_IO)  THEN
1504              CALL allgather_1i (LOCAL_CNT,ALL_PART_CNT,IERR)
1505     
1506              IF (myPE == PE_IO.AND.GLOBAL_CNT>0) THEN
1507                 WRITE(PVTU_UNIT,100) '      </PPointData>'
1508     
1509                 DO N = 0,NumPEs-1
1510                    IF(ALL_PART_CNT(N)>0) THEN
1511                       IF(TIME_DEPENDENT_FILENAME.AND.MODE==0) THEN
1512                          WRITE(VTU_NAME,20) TRIM(VTK_FILEBASE(VTK_REGION)),FRAME(VTK_REGION),N
1513                       ELSE
1514                          WRITE(VTU_NAME,25) TRIM(VTK_FILEBASE(VTK_REGION)),N
1515                       ENDIF
1516     
1517                       WRITE(PVTU_UNIT,110) '      <Piece Source="',TRIM(VTU_NAME),'"/>'
1518                    ENDIF
1519                 ENDDO
1520     
1521     
1522                 WRITE(PVTU_UNIT,100) '  </PPolyData>'
1523                 WRITE(PVTU_UNIT,100) '</VTKFile>'
1524                 CLOSE(PVTU_UNIT)
1525              ENDIF
1526           ENDIF
1527     
1528     20    FORMAT(A,"_",I4.4,"_",I5.5,".vtp")
1529     25    FORMAT(A,"_",I5.5,".vtp")
1530     
1531     100   FORMAT(A)
1532     110   FORMAT(A,A,A)
1533     
1534           RETURN
1535     
1536           END SUBROUTINE CLOSE_VTP_FILE_BIN
1537     
1538     
1539     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1540     !                                                                      C
1541     !  Module name: SETUP_VTK_REGION_PARTICLES                             C
1542     !                                                                      C
1543     !  Purpose: Filter the particles  based on the VTK region bounds and   C
1544     !           set the flag BELONGS_TO_VTK_SUBDOMAIN to .TRUE.            C
1545     !           to keep the particle.                                      C
1546     !                                                                      C
1547     !  Author: Jeff Dietiker                              Date: 11-Feb-15  C
1548     !  Reviewer:                                          Date:            C
1549     !                                                                      C
1550     !  Revision Number #                                  Date: ##-###-##  C
1551     !  Author: #                                                           C
1552     !  Purpose: #                                                          C
1553     !                                                                      C
1554     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1555           SUBROUTINE SETUP_VTK_REGION_PARTICLES
1556     
1557           USE vtk, only: VTK_REGION
1558           USE vtk, only: VTK_X_E, VTK_X_W, VTK_Y_S, VTK_Y_N, VTK_Z_B, VTK_Z_T
1559           USE vtk, only: VTK_NXS, VTK_NYS, VTK_NZS
1560           USE vtk, only: VTK_SLICE_TOL, VTK_SELECT_MODE
1561           USE vtk, only: BELONGS_TO_VTK_SUBDOMAIN
1562           USE discretelement, only: MAX_PIP,PIP,DES_POS_NEW
1563     
1564           IMPLICIT NONE
1565     
1566           INTEGER :: PC,LC1
1567           INTEGER :: NXS,NYS,NZS,NS
1568           INTEGER :: X_SLICE(DIM_I),Y_SLICE(DIM_J),Z_SLICE(DIM_K)
1569           DOUBLE PRECISION :: XE,XW,YS,YN,ZB,ZT
1570           DOUBLE PRECISION :: XP,YP,ZP,XP1,YP1,ZP1,XP2,YP2,ZP2,R
1571     
1572           DOUBLE PRECISION :: SLICE_TOL
1573           LOGICAL :: KEEP_XDIR,KEEP_YDIR,KEEP_ZDIR
1574     
1575     ! Variables related to gather
1576           integer lgathercnts(0:numpes-1), lproc
1577     
1578           CHARACTER(LEN=1) :: SELECT_PARTICLE_BY
1579     
1580     ! Get VTK region bounds
1581           XE = VTK_X_E(VTK_REGION)
1582           XW = VTK_X_W(VTK_REGION)
1583           YS = VTK_Y_S(VTK_REGION)
1584           YN = VTK_Y_N(VTK_REGION)
1585           ZB = VTK_Z_B(VTK_REGION)
1586           ZT = VTK_Z_T(VTK_REGION)
1587     
1588           NXS = VTK_NXS(VTK_REGION)
1589           NYS = VTK_NYS(VTK_REGION)
1590           NZS = VTK_NZS(VTK_REGION)
1591     
1592           SLICE_TOL = VTK_SLICE_TOL(VTK_REGION)
1593     
1594           SELECT_PARTICLE_BY = VTK_SELECT_MODE(VTK_REGION)
1595     
1596     ! get slice(s) location
1597           DO NS = 1,NXS
1598              X_SLICE(NS) = XW + (XE-XW)/(NXS-1)*(NS-1)
1599           ENDDO
1600     
1601           DO NS = 1,NYS
1602              Y_SLICE(NS) = YS + (YN-YS)/(NYS-1)*(NS-1)
1603           ENDDO
1604     
1605           DO NS = 1,NZS
1606              Z_SLICE(NS) = ZB + (ZT-ZB)/(NZS-1)*(NS-1)
1607           ENDDO
1608     
1609     
1610     
1611     ! Loop through all particles on local rank and keep a list of particles
1612     ! belonging to VTK region
1613     
1614           IF(ALLOCATED(BELONGS_TO_VTK_SUBDOMAIN)) DEALLOCATE(BELONGS_TO_VTK_SUBDOMAIN)
1615           ALLOCATE(BELONGS_TO_VTK_SUBDOMAIN(MAX_PIP))
1616     
1617           BELONGS_TO_VTK_SUBDOMAIN = .FALSE.
1618     
1619           LOCAL_CNT = 0
1620           PC = 1
1621           DO LC1 = 1, MAX_PIP
1622              IF(PC > PIP) EXIT
1623              IF(IS_NONEXISTENT(LC1)) CYCLE
1624              PC = PC+1
1625              IF(IS_GHOST(LC1) .OR. IS_ENTERING_GHOST(LC1) .OR. IS_EXITING_GHOST(LC1)) CYCLE
1626     
1627              SELECT CASE(SELECT_PARTICLE_BY)
1628                 CASE('C')  ! Particle center must be inside vtk region
1629     
1630                    XP = DES_POS_NEW(LC1,1)
1631                    YP = DES_POS_NEW(LC1,2)
1632                    ZP = DES_POS_NEW(LC1,3)
1633     
1634     ! X-direction
1635                    KEEP_XDIR=.FALSE.
1636                    IF(NXS==0) THEN
1637                       IF(XW<=XP.AND.XP<=XE) KEEP_XDIR=.TRUE.
1638                    ELSE
1639                       DO NS = 1,NXS
1640                          IF((X_SLICE(NS)-SLICE_TOL)<=XP.AND.XP<=(X_SLICE(NS)+SLICE_TOL)) KEEP_XDIR=.TRUE.
1641                       ENDDO
1642                    ENDIF
1643     
1644     ! Y-direction
1645                    KEEP_YDIR=.FALSE.
1646                    IF(NYS==0) THEN
1647                       IF(YS<=YP.AND.YP<=YN) KEEP_YDIR=.TRUE.
1648                    ELSE
1649                       DO NS = 1,NYS
1650                          IF((Y_SLICE(NS)-SLICE_TOL)<=YP.AND.YP<=(Y_SLICE(NS)+SLICE_TOL)) KEEP_YDIR=.TRUE.
1651                       ENDDO
1652                    ENDIF
1653     
1654     ! Z-direction
1655                    KEEP_ZDIR=.FALSE.
1656                    IF(NZS==0) THEN
1657                       IF(ZB<=ZP.AND.ZP<=ZT) KEEP_ZDIR=.TRUE.
1658                    ELSE
1659                       DO NS = 1,NZS
1660                          IF((Z_SLICE(NS)-SLICE_TOL)<=ZP.AND.ZP<=(Z_SLICE(NS)+SLICE_TOL)) KEEP_ZDIR=.TRUE.
1661                       ENDDO
1662                    ENDIF
1663     
1664     
1665                 CASE('P')  ! Entire particle must be inside vtk region
1666     
1667                    R = DES_RADIUS(LC1)
1668     
1669                    XP1 = DES_POS_NEW(LC1,1) - R
1670                    YP1 = DES_POS_NEW(LC1,2) - R
1671                    ZP1 = DES_POS_NEW(LC1,3) - R
1672     
1673                    XP2 = DES_POS_NEW(LC1,1) + R
1674                    YP2 = DES_POS_NEW(LC1,2) + R
1675                    ZP2 = DES_POS_NEW(LC1,3) + R
1676     
1677     ! X-direction
1678                    KEEP_XDIR=.FALSE.
1679                    IF(NXS==0) THEN
1680                       IF(XW<=XP1.AND.XP2<=XE) KEEP_XDIR=.TRUE.
1681                    ELSE
1682                       DO NS = 1,NXS
1683                          IF((X_SLICE(NS)-SLICE_TOL)<=XP1.AND.XP2<=(X_SLICE(NS)+SLICE_TOL)) KEEP_XDIR=.TRUE.
1684                       ENDDO
1685                    ENDIF
1686     
1687     ! Y-direction
1688                    KEEP_YDIR=.FALSE.
1689                    IF(NYS==0) THEN
1690                       IF(YS<=YP1.AND.YP2<=YN) KEEP_YDIR=.TRUE.
1691                    ELSE
1692                       DO NS = 1,NYS
1693                          IF((Y_SLICE(NS)-SLICE_TOL)<=YP1.AND.YP2<=(Y_SLICE(NS)+SLICE_TOL)) KEEP_YDIR=.TRUE.
1694                       ENDDO
1695                    ENDIF
1696     
1697     ! Z-direction
1698                    KEEP_ZDIR=.FALSE.
1699                    IF(NZS==0) THEN
1700                       IF(ZB<=ZP1.AND.ZP2<=ZT) KEEP_ZDIR=.TRUE.
1701                    ELSE
1702                       DO NS = 1,NZS
1703                          IF((Z_SLICE(NS)-SLICE_TOL)<=ZP1.AND.ZP2<=(Z_SLICE(NS)+SLICE_TOL)) KEEP_ZDIR=.TRUE.
1704                       ENDDO
1705                    ENDIF
1706     
1707     
1708                 CASE('I')  ! Particle must be inside or intersect the edge of the vtk region
1709     
1710                    R = DES_RADIUS(LC1)
1711     
1712                    XP1 = DES_POS_NEW(LC1,1) - R
1713                    YP1 = DES_POS_NEW(LC1,2) - R
1714                    ZP1 = DES_POS_NEW(LC1,3) - R
1715     
1716                    XP2 = DES_POS_NEW(LC1,1) + R
1717                    YP2 = DES_POS_NEW(LC1,2) + R
1718                    ZP2 = DES_POS_NEW(LC1,3) + R
1719     
1720     ! X-direction
1721                    KEEP_XDIR=.FALSE.
1722                    IF(NXS==0) THEN
1723                       IF(.NOT.(XE<=XP1.OR.XP2<=XW)) KEEP_XDIR=.TRUE.
1724                    ELSE
1725                       DO NS = 1,NXS
1726                          IF(.NOT.((X_SLICE(NS)+SLICE_TOL)<=XP1.OR.XP2<=(X_SLICE(NS)-SLICE_TOL))) KEEP_XDIR=.TRUE.
1727                       ENDDO
1728                    ENDIF
1729     
1730     ! Y-direction
1731                    KEEP_YDIR=.FALSE.
1732                    IF(NYS==0) THEN
1733                       IF(.NOT.(YN<=YP1.OR.YP2<=YS)) KEEP_YDIR=.TRUE.
1734                    ELSE
1735                       DO NS = 1,NYS
1736                          IF(.NOT.((Y_SLICE(NS)+SLICE_TOL)<=YP1.OR.YP2<=(Y_SLICE(NS)-SLICE_TOL))) KEEP_YDIR=.TRUE.
1737                       ENDDO
1738                    ENDIF
1739     
1740     ! Z-direction
1741                    KEEP_ZDIR=.FALSE.
1742                    IF(NZS==0) THEN
1743                       IF(.NOT.(ZT<=ZP1.OR.ZP2<=ZB)) KEEP_ZDIR=.TRUE.
1744                    ELSE
1745                       DO NS = 1,NZS
1746                          IF(.NOT.((Z_SLICE(NS)+SLICE_TOL)<=ZP1.OR.ZP2<=(Z_SLICE(NS)-SLICE_TOL))) KEEP_ZDIR=.TRUE.
1747                       ENDDO
1748                    ENDIF
1749     
1750     
1751                 CASE DEFAULT
1752                       print*,'should not be here'
1753              END SELECT
1754     
1755     ! Now combine
1756              IF(KEEP_XDIR.AND.KEEP_YDIR.AND.KEEP_ZDIR) THEN
1757                 BELONGS_TO_VTK_SUBDOMAIN(LC1) = .TRUE.
1758                 LOCAL_CNT = LOCAL_CNT + 1
1759              ENDIF
1760           ENDDO ! particle loop
1761     
1762     
1763     ! Calculate the total number of particles system-wide.
1764           call global_sum(LOCAL_CNT, GLOBAL_CNT)
1765     
1766     ! No need to set the send/reccv when using distributed IO
1767           IF (BDIST_IO) RETURN
1768     ! Set the send count from the local process.
1769           igath_sendcnt = LOCAL_CNT
1770     
1771     ! Collect the number of particles on each rank.all ranks.
1772           lgathercnts = 0
1773           lgathercnts(myPE) = LOCAL_CNT
1774           call global_sum(lgathercnts,igathercnts)
1775     
1776     ! Calculate the rank displacements.
1777           idispls(0) = 0
1778           DO lPROC = 1,NUMPEs-1
1779              idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
1780           ENDDO
1781     
1782           RETURN
1783     
1784           END SUBROUTINE SETUP_VTK_REGION_PARTICLES
1785     
1786           END MODULE VTP
1787