File: /nfs/home/0/users/jenkins/mfix.git/model/des/pic_bc_routines.f
1
2
3
4
5
6
7
8
9 SUBROUTINE PIC_APPLY_WALLBC_STL
10
11
12
13
14 use discretelement, only: DES_POS_NEW, DES_POS_OLD
15
16 use discretelement, only: DES_VEL_NEW
17
18 use discretelement, only: DES_RADIUS
19
20 use discretelement, only: PEA
21
22 use discretelement, only: MAX_PIP
23
24 use discretelement, only: cellneighbor_facet_num
25
26 use discretelement, only: cellneighbor_facet
27
28 use discretelement, only: DG_PIJK
29
30 use stl, only: NO_NEIGHBORING_FACET_DES
31
32 use stl, only: VERTEX
33
34 use stl, only: NORM_FACE
35
36 use pic_bc, only: PIC_BCMI, PIC_BCMO
37
38
39
40 use param1, only: ZERO, ONE, UNDEFINED
41
42 use discretelement, only: DIMN
43
44
45
46 use des_stl_functions, only: intersectLnPlane
47 use des_stl_functions, only: checkPTonTriangle
48
49 use error_manager
50
51 IMPLICIT NONE
52
53
54
55
56 INTEGER NF, LL
57
58 INTEGER :: CELL_COUNT
59
60 DOUBLE PRECISION, DIMENSION(DIMN) :: REF_LINE, DIR_LINE
61
62 DOUBLE PRECISION, DIMENSION(DIMN) :: REF_PLANE, NORM_PLANE
63
64 DOUBLE PRECISION :: LINE_T
65
66 LOGICAL :: ONTRIANGLE
67
68 DOUBLE PRECISION, DIMENSION(DIMN) :: POINT_ONPLANE
69
70
71
72 DOUBLE PRECISION :: dist_from_facet
73 DOUBLE PRECISION :: veldotnorm
74
75 CALL INIT_ERR_MSG("PIC_APPLY_WALLBC_STL")
76
77 DO LL = 1, MAX_PIP
78
79
80 IF(.NOT.PEA(LL,1)) CYCLE
81
82
83 IF (NO_NEIGHBORING_FACET_DES(DG_PIJK(LL))) CYCLE
84
85
86 CELLLOOP: DO CELL_COUNT=1, CELLNEIGHBOR_FACET_NUM(DG_PIJK(LL))
87
88
89 = CELLNEIGHBOR_FACET(DG_PIJK(LL))%P(CELL_COUNT)
90
91
92
93 = -UNDEFINED
94
95 ONTRIANGLE = .FALSE.
96
97
98
99
100 (1:DIMN) = DES_POS_OLD(:,LL)
101 DIR_LINE(1:DIMN) = DES_POS_NEW(:,LL) - DES_POS_OLD(:,LL)
102
103 NORM_PLANE = NORM_FACE(:,NF)
104 REF_PLANE = VERTEX(1,:,NF)
105
106 VELDOTNORM = DOT_PRODUCT(NORM_PLANE, DES_VEL_NEW(:,LL))
107
108
109
110 IF(VELDOTNORM.GT.0.d0) CYCLE
111
112 CALL INTERSECTLNPLANE(REF_LINE, DIR_LINE, REF_PLANE, &
113 NORM_PLANE, LINE_T)
114
115
116 IF(LINE_T.GE.ZERO.AND.LINE_T.LT.ONE) THEN
117
118 POINT_ONPLANE = REF_LINE + LINE_T*DIR_LINE
119
120
121
122
123 CALL CHECKPTONTRIANGLE(POINT_ONPLANE(:), VERTEX(:,:,NF),&
124 ONTRIANGLE)
125
126 IF(ONTRIANGLE) THEN
127
128 DIST_FROM_FACET = ABS(DOT_PRODUCT(DES_POS_NEW(:,LL) -&
129 POINT_ONPLANE(:), NORM_PLANE(:)))
130
131 DES_POS_NEW(:,LL) = POINT_ONPLANE + &
132 DIST_FROM_FACET*NORM_PLANE
133
134
135
136 CALL PIC_REFLECT_PART(LL, NORM_PLANE(:))
137
138 EXIT CELLLOOP
139 ENDIF
140 ENDIF
141
142 END DO CELLLOOP
143
144 END DO
145
146
147 IF(PIC_BCMI > 0) CALL PIC_MI_BC
148 IF(PIC_BCMO > 0) CALL PIC_MO_BC
149
150 CALL FINL_ERR_MSG
151
152 RETURN
153 END SUBROUTINE PIC_APPLY_WALLBC_STL
154
155
156
157
158
159
160
161
162
163
164
165
166 SUBROUTINE PIC_MO_BC
167
168 use discretelement
169 use pic_bc
170 use bc
171 USE error_manager
172 USE mpi_utility
173
174
175 implicit none
176
177 INTEGER :: IJK
178 INTEGER :: LC, LP, NP
179 INTEGER :: BCV, BCV_I
180
181 DOUBLE PRECISION :: DIST
182
183 INTEGER :: PIP_DEL_COUNT
184 INTEGER, dimension(0:numpes-1) :: del_count_all
185 INTEGER :: IPROC
186
187
188 CALL INIT_ERR_MSG("PIC_MO_BC")
189
190 PIP_DEL_COUNT = 0
191 DO BCV_I = 1, PIC_BCMO
192
193 BCV = PIC_BCMO_MAP(BCV_I)
194
195 SELECT CASE (BC_PLANE(BCV))
196 END SELECT
197
198 DO LC=PIC_BCMO_IJKSTART(BCV_I), PIC_BCMO_IJKEND(BCV_I)
199 IJK = PIC_BCMO_IJK(LC)
200 DO LP= 1,PINC(IJK)
201
202 NP = PIC(IJK)%p(LP)
203
204 IF(PEA(NP,4)) CYCLE
205
206 SELECT CASE (BC_PLANE(BCV))
207 CASE('S'); DIST = YN(BC_J_s(BCV)-1) - DES_POS_NEW(2,NP)
208 CASE('N'); DIST = DES_POS_NEW(2,NP) - YN(BC_J_s(BCV))
209 CASE('W'); DIST = XE(BC_I_w(BCV)-1) - DES_POS_NEW(1,NP)
210 CASE('E'); DIST = DES_POS_NEW(1,NP) - XE(BC_I_w(BCV))
211 CASE('B'); DIST = ZT(BC_K_b(BCV)-1) - DES_POS_NEW(3,NP)
212 CASE('T'); DIST = DES_POS_NEW(3,NP) - ZT(BC_K_b(BCV))
213 END SELECT
214
215 IF(DIST .lt. zero ) THEN
216
217
218
219 CALL DELETE_PARCEL(NP)
220 PIP_DEL_COUNT = PIP_DEL_COUNT+1
221 ENDIF
222
223 ENDDO
224 ENDDO
225 ENDDO
226
227
228 IF(PIC_REPORT_DELETION_STATS) then
229
230 del_count_all(:) = 0
231 del_count_all(mype) = pip_del_count
232 call global_all_sum(del_count_all(0:numpes-1))
233
234 IF(SUM(del_COUNT_ALL(:)).GT.0) THEN
235 WRITE(err_msg,'(/,2x,A,2x,i10)')&
236 'TOTAL NUMBER OF PARCELS DELETED GLOBALLY = ', SUM(DEL_COUNT_ALL(:))
237
238 call flush_err_msg(header = .false., footer = .false.)
239
240 DO IPROC = 0, NUMPES-1
241 IF(DMP_LOG) WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
242 'PARCELS ADDED ON PROC:', IPROC,' EQUAL TO', pip_del_count
243 ENDDO
244 ENDIF
245
246 ENDIF
247
248
249 CALL FINL_ERR_MSG
250 RETURN
251 END SUBROUTINE PIC_MO_BC
252
253
254
255
256
257
258
259
260
261
262
263
264 SUBROUTINE DELETE_PARCEL(NP)
265
266 USE compar
267 USE constant
268 USE des_bc
269 USE discretelement
270 USE funits
271 USE geometry
272 USE indices
273 USE param1
274 USE physprop
275 USE mfix_pic
276 USE functions
277
278 IMPLICIT NONE
279
280
281 INTEGER, INTENT(IN) :: NP
282
283 PEA(NP,:) = .FALSE.
284
285 DES_POS_OLD(:,NP) = ZERO
286 DES_POS_NEW(:,NP) = ZERO
287 DES_VEL_OLD(:,NP) = ZERO
288 DES_VEL_NEW(:,NP) = ZERO
289 OMEGA_OLD(:,NP) = ZERO
290 OMEGA_NEW(:,NP) = ZERO
291 DES_RADIUS(NP) = ZERO
292 PMASS(NP) = ZERO
293 PVOL(NP) = ZERO
294 RO_Sol(NP) = ZERO
295 OMOI(NP) = ZERO
296
297 DES_STAT_WT(NP) = ZERO
298
299 FC(:,NP) = ZERO
300
301 PIP = PIP - 1
302
303 RETURN
304 END SUBROUTINE DELETE_PARCEL
305
306
307
308
309
310
311
312
313 SUBROUTINE PIC_MI_BC
314
315 USE run
316 USE discretelement
317 USE geometry
318 USE constant
319 USE cutcell
320 USE mfix_pic
321 USE pic_bc
322 USE bc
323 USE mpi_utility
324 USE randomno
325 use error_manager
326 USE functions
327
328 IMPLICIT NONE
329
330
331
332 INTEGER :: WDIR, IDIM, IPCOUNT, LAST_EMPTY_SPOT, NEW_SPOT
333 INTEGER :: BCV, BCV_I, L, LC, PIP_ADD_COUNT, IPROC
334 INTEGER :: IFLUID, JFLUID, KFLUID, IJK, M
335 DOUBLE PRECISION :: CORD_START(3), DOML(3),WALL_NORM(3)
336 DOUBLE PRECISION :: AREA_INFLOW, VEL_INFLOW(DIMN), STAT_WT
337
338 LOGICAL :: DELETE_PART
339 DOUBLE PRECISION :: EPS_INFLOW(DES_MMAX)
340 DOUBLE PRECISION REAL_PARTS(DES_MMAX), COMP_PARTS(DES_MMAX), VEL_NORM_MAG, VOL_INFLOW, VOL_IJK
341
342
343
344 LOGICAL :: CONST_NPC, CONST_STATWT
345
346 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RANDPOS
347 INTEGER :: CNP_CELL_COUNT
348 DOUBLE PRECISION :: RNP_CELL_COUNT
349 integer :: lglobal_id
350 integer, dimension(0:numpes-1) :: add_count_all
351
352 LOGICAL, parameter :: setDBG = .true.
353 LOGICAL :: dFlag
354
355 type :: ty_spotlist
356 integer spot
357 type(ty_spotlist),pointer :: next => NULL()
358 end type ty_spotlist
359
360 type(ty_spotlist),pointer :: &
361 cur_spotlist => NULL(), &
362 prev_spotlist => NULL(), &
363 temp_spotlist => NULL()
364
365
366
367 CALL INIT_ERR_MSG("PIC_MI_BC")
368
369 dFlag = (DMP_LOG .AND. setDBG)
370 PIP_ADD_COUNT = 0
371 LAST_EMPTY_SPOT = 0
372
373 ALLOCATE(cur_spotlist); cur_spotlist%spot = -1
374
375 DO BCV_I = 1, PIC_BCMI
376 BCV = PIC_BCMI_MAP(BCV_I)
377
378 WALL_NORM(1:3) = PIC_BCMI_NORMDIR(BCV_I,1:3)
379
380
381 = 0
382 DO IDIM = 1, DIMN
383 WDIR = WDIR + ABS(WALL_NORM(IDIM))*IDIM
384 end DO
385
386 DO LC=PIC_BCMI_IJKSTART(BCV_I), PIC_BCMI_IJKEND(BCV_I)
387 IJK = PIC_BCMI_IJK(LC)
388
389 IF(.NOT.FLUID_AT(IJK)) CYCLE
390
391 IFLUID = I_OF(IJK)
392 JFLUID = J_OF(IJK)
393 KFLUID = K_OF(IJK)
394
395 CORD_START(1) = XE(IFLUID) - PIC_BCMI_OFFSET (BCV_I,1)*DX(IFLUID)
396
397 CORD_START(2) = YN(JFLUID) - PIC_BCMI_OFFSET (BCV_I,2)*DY(JFLUID)
398
399
400 CORD_START(3) = merge(zero, ZT(KFLUID) - PIC_BCMI_OFFSET (BCV_I,3)*DZ(KFLUID), no_k)
401
402 DOML(1) = DX(IFLUID)
403 DOML(2) = DY(JFLUID)
404 DOML(3) = MERGE(DZ(1), DZ(KFLUID), NO_K)
405
406 AREA_INFLOW = DOML(1)*DOML(2)*DOML(3)/DOML(WDIR)
407
408 VOL_IJK = DOML(1)*DOML(2)*DOML(3)
409
410 DOML(WDIR) = ZERO
411
412
413
414 DO M = 1, DES_MMAX
415
416 IF(SOLIDS_MODEL(M) /= 'PIC') CYCLE
417 EPS_INFLOW(M) = BC_ROP_S(BCV, M)/DES_RO_S(M)
418 VEL_INFLOW(1) = BC_U_S(BCV, M)
419 VEL_INFLOW(2) = BC_V_S(BCV, M)
420 VEL_INFLOW(3) = BC_W_S(BCV, M)
421
422 VEL_NORM_MAG = ABS(DOT_PRODUCT(VEL_INFLOW(1:DIMN), WALL_NORM(1:DIMN)))
423 VOL_INFLOW = AREA_INFLOW*VEL_NORM_MAG*DTSOLID
424
425 REAL_PARTS(M) = 6.d0*EPS_INFLOW(M)*VOL_INFLOW/(PI*(DES_D_p0(M)**3.d0))
426 COMP_PARTS(M) = zero
427
428 CONST_NPC = (BC_PIC_MI_CONST_NPC (BCV, M) .ne. 0)
429 CONST_STATWT = (BC_PIC_MI_CONST_STATWT(BCV, M) .ne. ZERO)
430 IF(CONST_NPC) THEN
431 IF(EPS_INFLOW(M).GT.ZERO) &
432 COMP_PARTS(M) = REAL(BC_PIC_MI_CONST_NPC(BCV, M))* &
433 VOL_INFLOW/VOL_IJK
434 ELSEIF(CONST_STATWT) THEN
435 COMP_PARTS(M) = REAL_PARTS(M)/ &
436 BC_PIC_MI_CONST_STATWT(BCV, M)
437 ENDIF
438
439 pic_bcmi_rnp(LC,M) = pic_bcmi_rnp(LC,M) + REAL_PARTS(M)
440
441 pic_bcmi_cnp(LC,M) = pic_bcmi_cnp(LC,M) + COMP_PARTS(M)
442
443
444 IF(pic_bcmi_cnp(LC,M).GE.1.d0) THEN
445 CNP_CELL_COUNT = INT(pic_bcmi_cnp(LC,M))
446 pic_bcmi_cnp(LC,M) = pic_bcmi_cnp(LC,M) - CNP_CELL_COUNT
447
448 RNP_CELL_COUNT = pic_bcmi_rnp(LC,M)
449 pic_bcmi_rnp(LC,M) = zero
450
451
452 = RNP_CELL_COUNT/REAL(CNP_CELL_COUNT)
453 ALLOCATE(RANDPOS(CNP_CELL_COUNT*DIMN))
454 CALL UNI_RNO(RANDPOS(:))
455
456 DO IPCOUNT = 1, CNP_CELL_COUNT
457
458
459 CALL PIC_FIND_EMPTY_SPOT(LAST_EMPTY_SPOT, NEW_SPOT)
460
461 DES_POS_OLD(1:DIMN, NEW_SPOT) = CORD_START(1:DIMN) &
462 & + RANDPOS((IPCOUNT-1)*DIMN+1: &
463 & (IPCOUNT-1)*DIMN+DIMN)*DOML(1:DIMN)
464 DES_POS_NEW(:, NEW_SPOT) = DES_POS_OLD(:,NEW_SPOT)
465 DES_VEL_OLD(1:DIMN,NEW_SPOT) = VEL_INFLOW(1:DIMN)
466
467 DES_VEL_NEW(:, NEW_SPOT) = DES_VEL_OLD(:, NEW_SPOT)
468
469 DES_RADIUS(NEW_SPOT) = DES_D_p0(M)*HALF
470
471 RO_Sol(NEW_SPOT) = DES_RO_S(M)
472
473 DES_STAT_WT(NEW_SPOT) = STAT_WT
474
475 PIJK(NEW_SPOT, 1) = IFLUID
476 PIJK(NEW_SPOT, 2) = JFLUID
477 PIJK(NEW_SPOT, 3) = KFLUID
478 PIJK(NEW_SPOT, 4) = IJK
479 PIJK(NEW_SPOT, 5) = M
480
481 PVOL(NEW_SPOT) = (4.0d0/3.0d0)*Pi*DES_RADIUS(NEW_SPOT)**3
482 PMASS(NEW_SPOT) = PVOL(NEW_SPOT)*RO_SOL(NEW_SPOT)
483
484 FC(:, NEW_SPOT) = zero
485 DELETE_PART = .false.
486 IF(PIC_BCMI_INCL_CUTCELL(BCV_I)) &
487 CALL CHECK_IF_PARCEL_OVELAPS_STL &
488 (des_pos_new(1:dimn, NEW_SPOT), &
489 DELETE_PART)
490
491 IF(.NOT.DELETE_PART) THEN
492
493 PIP = PIP+1
494 PIP_ADD_COUNT = PIP_ADD_COUNT + 1
495 PEA(NEW_SPOT, 1) = .true.
496 PEA(NEW_SPOT, 2:4) = .false.
497
498 ALLOCATE(temp_spotlist)
499 temp_spotlist%spot = new_spot
500 temp_spotlist%next => cur_spotlist
501 cur_spotlist => temp_spotlist
502 nullify(temp_spotlist)
503
504 ELSE
505 PEA(NEW_SPOT, 1) = .false.
506 LAST_EMPTY_SPOT = NEW_SPOT - 1
507 ENDIF
508
509
510
511
512
513
514
515
516 END DO
517 DEALLOCATE(RANDPOS)
518
519 end IF
520 end DO
521 end DO
522 end DO
523
524
525 (:) = 0
526 add_count_all(mype) = pip_add_count
527 call global_all_sum(add_count_all(0:numpes-1))
528 lglobal_id = imax_global_id + sum(add_count_all(0:mype-1))
529
530 do l = 1,pip_add_count
531 lglobal_id = lglobal_id + 1
532 iglobal_id(cur_spotlist%spot)= lglobal_id
533 prev_spotlist=> cur_spotlist
534 cur_spotlist => cur_spotlist%next
535 deallocate(prev_spotlist)
536 end do
537 deallocate(cur_spotlist)
538 imax_global_id = imax_global_id + sum(add_count_all(0:numpes-1))
539
540 IF(PIC_REPORT_SEEDING_STATS) then
541
542 IF(SUM(ADD_COUNT_ALL(:)).GT.0) THEN
543 WRITE(err_msg,'(/,2x,A,2x,i10)') &
544 'TOTAL NUMBER OF PARCELS ADDED GLOBALLY = ',&
545 SUM(ADD_COUNT_ALL(:))
546
547 call flush_err_msg(header = .false., footer = .false.)
548
549 DO IPROC = 0, NUMPES-1
550 IF(DMP_LOG) WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
551 'PARCELS ADDED ON PROC:', IPROC,&
552 ' EQUAL TO', ADD_COUNT_ALL(IPROC)
553 ENDDO
554 ENDIF
555
556 ENDIF
557
558 CALL FINL_ERR_MSG
559 RETURN
560 END SUBROUTINE PIC_MI_BC
561
562
563
564
565
566
567
568
569
570
571 SUBROUTINE PIC_FIND_EMPTY_SPOT(LAST_INDEX, EMPTY_SPOT)
572 USE funits
573 USE mpi_utility
574 USE error_manager
575 USE discretelement, only: max_pip, pea
576 IMPLICIT NONE
577
578
579
580 INTEGER, INTENT(INOUT) :: LAST_INDEX
581 INTEGER, INTENT(OUT) :: EMPTY_SPOT
582
583
584
585 LOGICAL :: SPOT_FOUND
586 INTEGER :: LL
587
588
589 CALL INIT_ERR_MSG("PIC_FIND_EMPTY_SPOT")
590
591 IF(LAST_INDEX.EQ.MAX_PIP) THEN
592 WRITE(ERR_MSG,2001)
593
594 CALL FLUSH_ERR_MSG(abort = .true.)
595 call mfix_exit(mype)
596 ENDIF
597 SPOT_FOUND = .false.
598
599 DO LL = LAST_INDEX+1, MAX_PIP
600
601 if(.NOT.PEA(LL,1)) THEN
602 EMPTY_SPOT = LL
603 LAST_INDEX = LL
604 SPOT_FOUND = .true.
605 EXIT
606 ENDIF
607 ENDDO
608
609 IF(.NOT.SPOT_FOUND) THEN
610 WRITE(ERR_MSG,2002)
611 CALL FLUSH_ERR_MSG(abort = .true.)
612 ENDIF
613
614 2001 FORMAT(/,5X, &
615 & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5X, &
616 & 'NO MORE EMPTY SPOT IN THE PARTICLE ARRAY TO ADD A NEW PARTICLE',/5X, &
617 & 'TERMINAL ERROR: STOPPING')
618
619 2002 FORMAT(/,5X, &
620 & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5X, &
621 & 'COULD NOT FIND A SPOT FOR ADDING NEW PARTICLE',/5X, &
622 & 'INCREASE THE SIZE OF THE INITIAL ARRAYS', 5X, &
623 & 'TERMINAL ERROR: STOPPING')
624
625 CALL FINL_ERR_MSG
626 END SUBROUTINE PIC_FIND_EMPTY_SPOT
627
628
629
630
631
632
633
634
635
636
637
638 SUBROUTINE PIC_REFLECT_PART(LL, WALL_NORM)
639
640 USE discretelement, only : dimn, DES_VEL_NEW
641 USE mfix_pic, only : MPPIC_COEFF_EN_WALL, &
642 & MPPIC_COEFF_ET_WALL
643 IMPLICIT NONE
644
645
646
647 INTEGER, INTENT(IN) :: LL
648 DOUBLE PRECISION, INTENT(IN) :: WALL_NORM(DIMN)
649
650
651
652
653 DOUBLE PRECISION :: VEL_NORMMAG_APP
654
655
656
657 DOUBLE PRECISION :: VEL_NORM_APP(DIMN), VEL_TANG_APP(DIMN)
658
659
660
661
662 DOUBLE PRECISION :: VEL_NORM_SEP(DIMN), VEL_TANG_SEP(DIMN)
663
664 DOUBLE PRECISION :: COEFF_REST_EN, COEFF_REST_ET
665
666
667
668 = DOT_PRODUCT(WALL_NORM(1:DIMN), DES_VEL_NEW(:, LL))
669
670
671
672
673 (1:DIMN) = VEL_NORMMAG_APP*WALL_NORM(1:DIMN)
674
675 VEL_TANG_APP(:) = DES_VEL_NEW(:, LL) - VEL_NORM_APP(1:DIMN)
676
677
678
679
680
681
682 = MPPIC_COEFF_EN_WALL
683 COEFF_REST_ET = MPPIC_COEFF_ET_WALL
684
685
686 (1:DIMN) = -COEFF_REST_EN*VEL_NORM_APP(1:DIMN)
687
688 VEL_TANG_SEP(1:DIMN) = COEFF_REST_ET*VEL_TANG_APP(1:DIMN)
689
690 DES_VEL_NEW(:, LL) = VEL_NORM_SEP(1:DIMN) + VEL_TANG_SEP(1:DIMN)
691
692 END SUBROUTINE PIC_REFLECT_PART
693
694
695
696
697
698
699
700
701
702
703
704 SUBROUTINE PIC_FIND_NEW_CELL(LL)
705 USE discretelement
706 use mpi_utility
707 USE cutcell
708 USE functions
709
710 IMPLICIT NONE
711
712
713
714
715 INTEGER, INTENT(IN) :: LL
716
717
718
719 DOUBLE PRECISION :: XPOS, YPOS, ZPOS
720 INTEGER :: I, J, K
721
722
723 = PIJK(LL,1)
724 J = PIJK(LL,2)
725 K = PIJK(LL,3)
726 XPOS = DES_POS_NEW(1,LL)
727 YPOS = DES_POS_NEW(2,LL)
728 IF (DIMN .EQ. 3) THEN
729 ZPOS = DES_POS_NEW(3,LL)
730 ENDIF
731
732 IF(XPOS >= XE(I-1) .AND. XPOS < XE(I)) THEN
733 PIJK(LL,1) = I
734 ELSEIF(XPOS >= XE(I)) THEN
735 PIJK(LL,1) = I+1
736 ELSE
737 PIJK(LL,1) = I-1
738 END IF
739
740 IF(YPOS >= YN(J-1) .AND. YPOS < YN(J)) THEN
741 PIJK(LL,2) = J
742 ELSEIF(YPOS >= YN(J))THEN
743 PIJK(LL,2) = J+1
744 ELSE
745 PIJK(LL,2) = J-1
746 END IF
747
748 IF(DIMN.EQ.2) THEN
749 PIJK(LL,3) = 1
750 ELSE
751 IF(ZPOS >= ZT(K-1) .AND. ZPOS < ZT(K)) THEN
752 PIJK(LL,3) = K
753 ELSEIF(ZPOS >= ZT(K)) THEN
754 PIJK(LL,3) = K+1
755 ELSE
756 PIJK(LL,3) = K-1
757 END IF
758 ENDIF
759
760 I = PIJK(LL,1)
761 J = PIJK(LL,2)
762 K = PIJK(LL,3)
763
764 IF(I.EQ.IEND1+1) then
765 IF(XPOS >= XE(IEND1-1) .AND. XPOS <= XE(IEND1)) PIJK(LL,1) = IEND1
766 ENDIF
767
768 IF(J.EQ.JEND1+1) then
769 IF(YPOS >= YN(JEND1-1) .AND. YPOS <= YN(JEND1)) PIJK(LL,2) = JEND1
770 ENDIF
771
772 IF(DIMN.EQ.3.AND.K.EQ.KEND1+1) THEN
773 IF(ZPOS >= ZT(KEND1-1) .AND. ZPOS <= ZT(KEND1)) PIJK(LL,3) = KEND1
774 ENDIF
775
776 PIJK(LL,4) = FUNIJK(PIJK(LL,1), PIJK(LL,2),PIJK(LL,3))
777
778 END SUBROUTINE PIC_FIND_NEW_CELL
779
780
781
782
783
784
785
786
787
788
789
790
791
792 SUBROUTINE CHECK_IF_PARCEL_OVELAPS_STL(POSITION, &
793 OVERLAP_EXISTS)
794
795 USE constant
796 USE cutcell
797 USE des_stl_functions
798 USE desgrid
799 USE discretelement, only: dimn
800 USE functions
801 USE geometry
802 USE indices
803 USE mpi_utility
804 USE param1
805 USE run
806 USE stl
807 Implicit none
808
809 DOUBLE PRECISION, INTENT(IN) :: POSITION(DIMN)
810 LOGICAL, INTENT(OUT) :: OVERLAP_EXISTS
811
812 INTEGER I, J, K, IJK, NF, FOCUS_PARTICLE
813
814 INTEGER :: COUNT_FAC, COUNT, NEIGH_CELLS, &
815 NEIGH_CELLS_NONNAT, &
816 LIST_OF_CELLS(27), CELL_ID, I_CELL, J_CELL, K_CELL, cell_count , &
817 IMINUS1, IPLUS1, JMINUS1, JPLUS1, KMINUS1, KPLUS1
818
819 double precision :: velocity(dimn)
820
821 double precision, dimension(dimn) :: ref_line, dir_line
822
823 double precision, dimension(dimn) :: ref_plane, norm_plane
824
825 double precision :: line_t
826
827 logical :: ontriangle
828
829 double precision, dimension(dimn) :: point_onplane
830
831 FOCUS_PARTICLE = -1
832
833 OVERLAP_EXISTS = .false.
834
835 I_CELL = MIN(DG_IEND2,MAX(DG_ISTART2,IOFPOS(POSITION(1))))
836 J_CELL = MIN(DG_JEND2,MAX(DG_JSTART2,JOFPOS(POSITION(2))))
837 K_CELL = 1
838 IF(DO_K) K_CELL =MIN(DG_KEND2,MAX(DG_KSTART2,KOFPOS(POSITION(3))))
839
840 CELL_ID = DG_FUNIJK(I_CELL, J_CELL, K_CELL)
841 IF (NO_NEIGHBORING_FACET_DES(CELL_ID)) RETURN
842
843 LIST_OF_CELLS(:) = -1
844 NEIGH_CELLS = 0
845 NEIGH_CELLS_NONNAT = 0
846
847 COUNT_FAC = LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
848
849 IF (COUNT_FAC.gt.0) then
850
851 = NEIGH_CELLS + 1
852 LIST_OF_CELLS(NEIGH_CELLS) = CELL_ID
853 ENDIF
854
855 IPLUS1 = MIN( I_CELL + 1, DG_IEND2)
856 IMINUS1 = MAX( I_CELL - 1, DG_ISTART2)
857
858 JPLUS1 = MIN (J_CELL + 1, DG_JEND2)
859 JMINUS1 = MAX( J_CELL - 1, DG_JSTART2)
860
861 KPLUS1 = MIN (K_CELL + 1, DG_KEND2)
862 KMINUS1 = MAX( K_CELL - 1, DG_KSTART2)
863
864 DO K = KMINUS1, KPLUS1
865 DO J = JMINUS1, JPLUS1
866 DO I = IMINUS1, IPLUS1
867 IJK = DG_FUNIJK(I,J,K)
868 COUNT_FAC = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
869 IF(COUNT_FAC.EQ.0) CYCLE
870 NEIGH_CELLS_NONNAT = NEIGH_CELLS_NONNAT + 1
871 NEIGH_CELLS = NEIGH_CELLS + 1
872 LIST_OF_CELLS(NEIGH_CELLS) = IJK
873
874 ENDDO
875 ENDDO
876 ENDDO
877
878
879 DO CELL_COUNT = 1, NEIGH_CELLS
880 IJK = LIST_OF_CELLS(CELL_COUNT)
881
882 DO COUNT = 1, LIST_FACET_AT_DES(IJK)%COUNT_FACETS
883 NF = LIST_FACET_AT_DES(IJK)%FACET_LIST(COUNT)
884 line_t = -Undefined
885
886
887 = .false.
888
889
890
891
892
893
894 (1:dimn) = position(1:dimn)
895 dir_line(1:dimn) = NORM_FACE(1:dimn,NF)
896
897
898
899
900
901 (1:dimn) = NORM_FACE(1:dimn,NF)
902 ref_plane(1:dimn) = VERTEX(1, 1:dimn,NF)
903 CALL intersectLnPlane(ref_line, dir_line, ref_plane, &
904 norm_plane, line_t)
905 if(line_t.gt.zero) then
906
907
908
909 (1:dimn) = ref_line(1:dimn) + &
910 line_t*dir_line(1:dimn)
911
912
913
914
915 CALL checkPTonTriangle(point_onplane(1:dimn), &
916 VERTEX(:,:,NF), ontriangle)
917
918 if(ontriangle) then
919 OVERLAP_EXISTS = .true.
920
921
922
923 RETURN
924 endif
925 endif
926
927 ENDDO
928
929 end DO
930
931 RETURN
932
933 END SUBROUTINE CHECK_IF_PARCEL_OVELAPS_STL
934
935
936 SUBROUTINE write_this_facet_and_parcel(FID, position, velocity)
937 USE run
938 USE param1
939 USE discretelement
940 USE geometry
941 USE compar
942 USE constant
943 USE cutcell
944 USE funits
945 USE indices
946 USE physprop
947 USE parallel
948 USE stl
949 USE des_stl_functions
950 Implicit none
951
952 double precision, intent(in), dimension(dimn) :: position, velocity
953 Integer, intent(in) :: fid
954 Integer :: stl_unit, vtp_unit , k
955 CHARACTER(LEN=100) :: stl_fname, vtp_fname
956 real :: temp_array(3)
957
958 stl_unit = 1001
959 vtp_unit = 1002
960
961 WRITE(vtp_fname,'(A,"_OFFENDING_PARTICLE",".vtp")') TRIM(RUN_NAME)
962 WRITE(stl_fname,'(A,"_STL_FACE",".stl")') TRIM(RUN_NAME)
963
964 open(vtp_unit, file = vtp_fname, form='formatted')
965 open(stl_unit, file = stl_fname, form='formatted')
966
967 write(vtp_unit,"(a)") '<?xml version="1.0"?>'
968 write(vtp_unit,"(a,es24.16,a)") '<!-- time =',s_time,'s -->'
969 write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
970 ' version="0.1" byte_order="LittleEndian">'
971 write(vtp_unit,"(3x,a)") '<PolyData>'
972 write(vtp_unit,"(6x,a,i10.10,a,a)")&
973 '<Piece NumberOfPoints="',1,'" NumberOfVerts="0" ',&
974 'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
975 write(vtp_unit,"(9x,a)")&
976 '<PointData Scalars="Diameter" Vectors="Velocity">'
977 write(vtp_unit,"(12x,a)")&
978 '<DataArray type="Float32" Name="Diameter" format="ascii">'
979 write (vtp_unit,"(15x,es13.6)") (1.d0)
980 write(vtp_unit,"(12x,a)") '</DataArray>'
981
982 temp_array = zero
983 temp_array(1:DIMN) = velocity(1:dimn)
984 write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
985 'Name="Velocity" NumberOfComponents="3" format="ascii">'
986 write (vtp_unit,"(15x,3(es13.6,3x))")&
987 ((temp_array(k)),k=1,3)
988 write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
989
990 write(vtp_unit,"(9x,a)") '<CellData></CellData>'
991
992 temp_array = zero
993 temp_array(1:dimn) = position(1:dimn)
994 write(vtp_unit,"(9x,a)") '<Points>'
995 write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
996 'Name="Position" NumberOfComponents="3" format="ascii">'
997 write (vtp_unit,"(15x,3(es13.6,3x))")&
998 ((temp_array(k)),k=1,3)
999 write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
1000
1001 write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
1002 '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
1003 write(vtp_unit,"(6x,a,/3x,a,/a)")&
1004 '</Piece>','</PolyData>','</VTKFile>'
1005
1006
1007
1008 write(stl_unit,*)'solid vcg'
1009
1010 write(stl_unit,*) ' facet normal ', NORM_FACE(1:3,FID)
1011 write(stl_unit,*) ' outer loop'
1012 write(stl_unit,*) ' vertex ', VERTEX(1,1:3,FID)
1013 write(stl_unit,*) ' vertex ', VERTEX(2,1:3,FID)
1014 write(stl_unit,*) ' vertex ', VERTEX(3,1:3,FID)
1015 write(stl_unit,*) ' endloop'
1016 write(stl_unit,*) ' endfacet'
1017
1018 write(stl_unit,*)'endsolid vcg'
1019
1020 close(vtp_unit, status = 'keep')
1021 close(stl_unit, status = 'keep')
1022 write(*,*) 'wrote a facet and a parcel. now waiting'
1023 read(*,*)
1024 end SUBROUTINE write_this_facet_and_parcel
1025
1026