File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/cut_cell_preprocessing.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 SUBROUTINE CUT_CELL_PREPROCESSING
19
20 USE param
21 USE param1
22 USE parallel
23 USE constant
24 USE run
25 USE toleranc
26 USE geometry
27 USE indices
28 USE compar
29 USE sendrecv
30 USE quadric
31 USE cutcell
32 USE vtk
33 use cdist
34 USE fldvar
35 USE functions
36
37 IMPLICIT NONE
38
39 INTEGER :: SAFE_MODE_COUNT
40 DOUBLE PRECISION :: CPU_PP_START,CPU_PP_END
41
42 IF(.NOT.CG_HEADER_WAS_PRINTED) CALL PRINT_CG_HEADER
43
44 CALL CPU_TIME (CPU_PP_START)
45
46 CALL OPEN_CUT_CELL_FILES
47
48 CALL ALLOCATE_CUT_CELL_ARRAYS
49
50 CALL DEFINE_QUADRICS
51
52 CALL SET_3D_CUT_CELL_FLAGS
53
54 CALL GATHER_DATA
55
56
57
58
59
60 IF(WRITE_VTK_FILES.AND.(.NOT.BDIST_IO)) THEN
61 CALL WRITE_CUT_SURFACE_VTK
62 ENDIF
63
64
65
66 CALL SET_3D_CUT_U_CELL_FLAGS
67 CALL SET_3D_CUT_V_CELL_FLAGS
68 IF(DO_K) CALL SET_3D_CUT_W_CELL_FLAGS
69
70 CALL SET_3D_CUT_CELL_TREATMENT_FLAGS
71
72 CALL GET_3D_ALPHA_U_CUT_CELL
73 CALL GET_3D_ALPHA_V_CUT_CELL
74 IF(DO_K) CALL GET_3D_ALPHA_W_CUT_CELL
75
76 CALL SET_GHOST_CELL_FLAGS
77
78 CALL SET_ODXYZ_U_CUT_CELL
79 CALL SET_ODXYZ_V_CUT_CELL
80 IF(DO_K) CALL SET_ODXYZ_W_CUT_CELL
81
82 CALL GET_U_MASTER_CELLS
83 CALL GET_V_MASTER_CELLS
84 IF(DO_K) CALL GET_W_MASTER_CELLS
85
86 CALL SEND_RECEIVE_CUT_CELL_VARIABLES
87
88 CALL GET_DISTANCE_TO_WALL
89
90 CALL PRINT_GRID_STATISTICS
91
92 CALL CG_GET_BC_AREA
93
94 CALL FLOW_TO_VEL(.FALSE.)
95
96 CALL CG_FLOW_TO_VEL
97
98 CALL CONVERT_CG_MI_TO_PS
99
100
101 CALL CPU_TIME (CPU_PP_END)
102
103
104
105 IF(myPE == PE_IO) THEN
106 WRITE(*,20)'CARTESIAN GRID PRE-PROCESSING COMPLETED IN ',CPU_PP_END - CPU_PP_START, ' SECONDS.'
107 WRITE(*,10)'============================================================================='
108 ENDIF
109
110 IF(myPE == PE_IO) THEN
111
112 SAFE_MODE_COUNT = SUM(CG_SAFE_MODE)
113
114 IF(SAFE_MODE_COUNT>0) THEN
115
116
117 WRITE(*,10)'######################################################################'
118 WRITE(*,10)'######################################################################'
119 WRITE(*,10)'## ##'
120 WRITE(*,10)'## || ##'
121 WRITE(*,10)'## || ##'
122 WRITE(*,10)'## \/ ##'
123 WRITE(*,10)'## ##'
124 WRITE(*,10)'## ===> WARNING: RUNNING CARTESIAN GRID IN SAFE MODE ! <=== ##'
125 WRITE(*,10)'## ##'
126 WRITE(*,10)'## SAFE MODE ACTIVATED FOR : ##'
127 IF(CG_SAFE_MODE(1)==1) WRITE(*,10)'## - All scalar quantities ##'
128 IF(CG_SAFE_MODE(3)==1) WRITE(*,10)'## - X-Velocity (Gas and Solids) ##'
129 IF(CG_SAFE_MODE(4)==1) WRITE(*,10)'## - Y-Velocity (Gas and Solids) ##'
130 IF(CG_SAFE_MODE(5)==1) WRITE(*,10)'## - Z-Velocity (Gas and Solids) ##'
131 WRITE(*,10)'## ##'
132 WRITE(*,10)'## /\ ##'
133 WRITE(*,10)'## || ##'
134 WRITE(*,10)'## || ##'
135 WRITE(*,10)'## ##'
136 WRITE(*,10)'######################################################################'
137 WRITE(*,10)'######################################################################'
138
139
140
141 ENDIF
142 ENDIF
143
144
145 RETURN
146
147 10 FORMAT(1X,A)
148 20 FORMAT(1X,A,F8.2,A)
149
150 END SUBROUTINE CUT_CELL_PREPROCESSING
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167 SUBROUTINE PRINT_CG_HEADER
168
169 USE param
170 USE param1
171 USE parallel
172 USE compar
173 USE cutcell
174
175 IMPLICIT NONE
176
177
178
179
180 IF(myPE == PE_IO) THEN
181
182 WRITE(*,10)'============================================================================='
183 WRITE(*,10)' ____ ___________ _ ____ ____ __________ __________ '
184 WRITE(*,10)' | \ / | (_) \ \ / / | | | | '
185 WRITE(*,10)' | \ / ______| ___ \ \ / / | ______| | ______| '
186 WRITE(*,10)' | \ / |______ | | \ \/ / | | | | '
187 WRITE(*,10)' | \ / | | | \ / === | | | | ____ '
188 WRITE(*,10)' | |\ \/ /| ______| | | / \ | | | | |_ | '
189 WRITE(*,10)' | | \ / | | | | / /\ \ | |______ | |___| | '
190 WRITE(*,10)' | | \ / | | | | / / \ \ | | | | '
191 WRITE(*,10)' |___| \__/ |___| |___| /___/ \___\ |__________| |__________| '
192 WRITE(*,10)' '
193 WRITE(*,10)'============================================================================='
194 WRITE(*,10)'MFIX WITH CARTESIAN GRID IMPLEMENTATION.'
195
196
197 IF(RE_INDEXING) THEN
198 WRITE(*,10)'RE-INDEXING IS TURNED ON.'
199
200
201
202
203
204
205 ELSE
206 WRITE(*,10)'RE-INDEXING IS TURNED OFF.'
207 ENDIF
208
209
210
211 ENDIF
212
213 CG_HEADER_WAS_PRINTED = .TRUE.
214
215 10 FORMAT(1X,A)
216
217 RETURN
218
219 END SUBROUTINE PRINT_CG_HEADER
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246 SUBROUTINE EVAL_F(METHOD,x1,x2,x3,Q,f,CLIP_FLAG)
247
248 USE parallel
249 USE compar
250 USE sendrecv
251 USE quadric
252 USE cutcell
253
254 USE quadric
255
256
257 IMPLICIT NONE
258
259 DOUBLE PRECISION x1,x2,x3
260 DOUBLE PRECISION f
261 INTEGER :: Q,Q_ID,BCID
262 LOGICAL :: CLIP_FLAG
263 CHARACTER (LEN = 7) :: METHOD
264 CHARACTER(LEN=9) :: GR
265
266 INTEGER :: GROUP,GS,P
267
268 DOUBLE PRECISION,DIMENSION(DIM_GROUP,0:DIM_QUADRIC) :: F_G
269
270 SELECT CASE(METHOD)
271
272 CASE('QUADRIC')
273
274 F_G = - UNDEFINED
275
276 DO GROUP = 1, N_GROUP
277 GS = GROUP_SIZE(GROUP)
278 GR = TRIM(GROUP_RELATION(GROUP))
279
280 DO P = 1 , GS
281 Q_ID = GROUP_Q(GROUP,P)
282 CALL GET_F_QUADRIC(x1,x2,x3,Q_ID,F_G(GROUP,P),CLIP_FLAG)
283 ENDDO
284 IF(GR == 'AND') THEN
285 F_G(GROUP,0) = MAXVAL(F_G(GROUP,1:GS))
286 ELSEIF(GR == 'OR') THEN
287 F_G(GROUP,0) = MINVAL(F_G(GROUP,1:GS))
288 ELSEIF(GR == 'PIECEWISE') THEN
289 CALL REASSSIGN_QUADRIC(x1,x2,x3,GROUP,Q_ID)
290
291 CALL GET_F_QUADRIC(x1,x2,x3,Q_ID,F_G(GROUP,0),CLIP_FLAG)
292
293 ENDIF
294
295 ENDDO
296
297 f = F_G(1,0)
298
299 DO GROUP = 2, N_GROUP
300
301 GR = TRIM(RELATION_WITH_PREVIOUS(GROUP))
302
303 IF(GR =='AND') THEN
304 f = DMAX1(f,F_G(GROUP,0))
305 ELSEIF(GR =='OR') THEN
306 f = DMIN1(f,F_G(GROUP,0))
307 ENDIF
308
309 ENDDO
310
311
312 CASE('POLYGON')
313
314 CALL EVAL_POLY_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
315
316 CASE('USR_DEF')
317
318 CALL EVAL_USR_FCT(x1,x2,x3,Q,f,CLIP_FLAG)
319
320
321
322
323
324 CASE DEFAULT
325
326 WRITE(*,*)'ERROR IN SUBROUTINE EVAL_F.'
327 WRITE(*,*)'UNKNOWN METHOD:',METHOD
328 WRITE(*,*)'ACCEPTABLE METHODS:'
329 WRITE(*,*)'QUADRIC'
330 WRITE(*,*)'POLYGON'
331 WRITE(*,*)'USR_DEF'
332
333 CALL MFIX_EXIT(myPE)
334
335 END SELECT
336
337
338
339 RETURN
340 END SUBROUTINE EVAL_F
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358 SUBROUTINE INTERSECT_LINE(METHOD,xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
359
360 USE param
361 USE param1
362 USE parallel
363 USE constant
364 USE run
365 USE toleranc
366 USE geometry
367 USE indices
368 USE compar
369 USE sendrecv
370 USE quadric
371
372 IMPLICIT NONE
373 DOUBLE PRECISION:: x1,y1,z1,x2,y2,z2,x3,y3,z3
374 DOUBLE PRECISION:: xa,ya,za,xb,yb,zb,xc,yc,zc
375 INTEGER :: Q_ID,niter
376 DOUBLE PRECISION :: x_intersection
377 DOUBLE PRECISION :: f1,f2,f3,fa,fb
378 DOUBLE PRECISION :: t1,t2,t3
379 LOGICAL :: CLIP_FLAG,CLIP_FLAG1,CLIP_FLAG2,CLIP_FLAG3,INTERSECT_FLAG
380 CHARACTER (LEN=7) ::METHOD
381
382 x1 = xa
383 = ya
384 z1 = za
385 t1 = ZERO
386
387 x2 = xb
388 y2 = yb
389 z2 = zb
390 t2 = ONE
391
392 CALL EVAL_F(METHOD,x1,y1,z1,Q_ID,f1,CLIP_FLAG1)
393 CALL EVAL_F(METHOD,x2,y2,z2,Q_ID,f2,CLIP_FLAG2)
394
395
396
397
398
399
400 = 1
401
402 CLIP_FLAG = (CLIP_FLAG1).AND.(CLIP_FLAG2)
403
404 if(DABS(f1)<TOL_F) then
405 = UNDEFINED
406 yc = UNDEFINED
407 zc = UNDEFINED
408 INTERSECT_FLAG = .FALSE.
409 elseif(DABS(f2)<TOL_F) then
410 = UNDEFINED
411 yc = UNDEFINED
412 zc = UNDEFINED
413 INTERSECT_FLAG = .FALSE.
414 elseif(f1*f2 < ZERO) then
415 niter = 0
416 f3 = 2.0d0*TOL_F
417 do while ( (abs(f3) > TOL_F) .AND. (niter<ITERMAX_INT) )
418
419 t3 = t1 - f1*(t2-t1)/(f2-f1)
420
421 x3 = x1 + t3 * (x2 - x1)
422 y3 = y1 + t3 * (y2 - y1)
423 z3 = z1 + t3 * (z2 - z1)
424
425 CALL EVAL_F(METHOD,x3,y3,z3,Q_ID,f3,CLIP_FLAG3)
426 if(f1*f3<0) then
427 t2 = t3
428 f2 = f3
429 else
430 t1 = t3
431 f1 = f3
432 endif
433 niter = niter + 1
434
435 end do
436 if (niter < ITERMAX_INT) then
437 xc = x3
438 yc = y3
439 zc = z3
440 INTERSECT_FLAG = .TRUE.
441 else
442 WRITE(*,*)' Subroutine intersect_line:'
443 WRITE(*,*) 'Unable to find the intersection of quadric:',Q_ID
444 WRITE(*,1000)'between (x1,y1,z1)= ', xa,ya,za
445 WRITE(*,1000)' and (x2,y2,z2)= ', xb,yb,zb
446 CALL EVAL_F(METHOD,xa,ya,za,Q_ID,fa,CLIP_FLAG1)
447 CALL EVAL_F(METHOD,xb,yb,zb,Q_ID,fb,CLIP_FLAG1)
448 WRITE(*,1000)'f(x1,y1,z1) = ', fa
449 WRITE(*,1000)'f(x2,y2,z2) = ', fb
450 WRITE(*,1000)'Current Location (x3,y3,z3)= ', x3,y3,z3
451 WRITE(*,1000)'Current value of abs(f) = ', DABS(f3)
452 WRITE(*,1000)'Tolerance = ', TOL_F
453 WRITE(*,*)'Maximum number of iterations = ', ITERMAX_INT
454 WRITE(*,*) 'Please increase the intersection tolerance, '
455 WRITE(*,*) 'or the maximum number of iterations, and try again.'
456 WRITE(*,*) 'MFiX will exit now.'
457 CALL MFIX_EXIT(myPE)
458 x_intersection = UNDEFINED
459 INTERSECT_FLAG = .FALSE.
460
461 endif
462 else
463 xc = UNDEFINED
464 yc = UNDEFINED
465 zc = UNDEFINED
466 INTERSECT_FLAG = .FALSE.
467 endif
468
469 1000 FORMAT(A,3(2X,G12.5))
470
471
472 RETURN
473
474 END SUBROUTINE INTERSECT_LINE
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489 SUBROUTINE INTERSECT(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
490
491 USE param
492 USE param1
493 USE parallel
494 USE constant
495 USE run
496 USE toleranc
497 USE geometry
498 USE indices
499 USE compar
500 USE sendrecv
501 USE quadric
502 USE cutcell
503 USE polygon
504 USE functions
505
506 IMPLICIT NONE
507 CHARACTER (LEN=*) :: TYPE_OF_CELL
508 INTEGER :: IJK,I,J,K,Q_ID,N_int_x,N_int_y,N_int_z,N_USR
509 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb,xc,yc,zc
510 DOUBLE PRECISION :: Xi,Yi,Zi,Xc_backup,Yc_backup,Zc_backup
511 LOGICAL :: INTERSECT_FLAG
512
513 Xi = UNDEFINED
514 Yi = UNDEFINED
515 Zi = UNDEFINED
516
517
518
519
520
521 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
522
523
524
525
526 = X_NODE(7)
527 ya = Y_NODE(7)
528 za = Z_NODE(7)
529
530 xb = X_NODE(8)
531 yb = Y_NODE(8)
532 zb = Z_NODE(8)
533
534
535 N_int_x = 0
536 INTERSECT_X(IJK) = .FALSE.
537 Q_ID = 1
538 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
539 IF ((INTERSECT_FLAG).AND.(xc/=Xi)) THEN
540 N_int_x = N_int_x + 1
541 INTERSECT_X(IJK) = .TRUE.
542 xc_backup = Xi
543 Xi = xc
544 ENDIF
545
546 IF(N_int_x /= 1) THEN
547 Xi = UNDEFINED
548 INTERSECT_X(IJK) = .FALSE.
549 ENDIF
550
551
552 DO Q_ID = 1, N_POLYGON
553 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_X(IJK),xc,yc,zc)
554 IF(INTERSECT_X(IJK)) Xi = xc
555 ENDDO
556
557 DO N_USR= 1, N_USR_DEF
558 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_X(IJK),xc,yc,zc)
559 IF(INTERSECT_X(IJK)) Xi = xc
560 ENDDO
561
562
563
564
565
566
567 IF(TYPE_OF_CELL=='U_MOMENTUM') THEN
568 IF(SNAP(IJK)) THEN
569 INTERSECT_X(IJK) = .TRUE.
570 I = I_OF(IJK)
571 Xi = XG_E(I)
572 ENDIF
573 ENDIF
574
575
576
577
578
579
580 = X_NODE(6)
581 ya = Y_NODE(6)
582 za = Z_NODE(6)
583
584 N_int_y = 0
585 INTERSECT_Y(IJK) = .FALSE.
586 Q_ID = 1
587 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
588 IF ((INTERSECT_FLAG).AND.(yc/=Yi)) THEN
589 N_int_y = N_int_y + 1
590 INTERSECT_Y(IJK) = .TRUE.
591 yc_backup = Yi
592 Yi = yc
593 ENDIF
594
595 IF(N_int_y /= 1) THEN
596 Yi = UNDEFINED
597 INTERSECT_Y(IJK) = .FALSE.
598 ENDIF
599
600 DO Q_ID = 1, N_POLYGON
601 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Y(IJK),xc,yc,zc)
602 IF(INTERSECT_Y(IJK)) Yi = yc
603 ENDDO
604
605 DO N_USR= 1, N_USR_DEF
606 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Y(IJK),xc,yc,zc)
607 IF(INTERSECT_Y(IJK)) Yi = yc
608 ENDDO
609
610
611
612
613
614
615 IF(TYPE_OF_CELL=='V_MOMENTUM') THEN
616 IF(SNAP(IJK)) THEN
617 INTERSECT_Y(IJK) = .TRUE.
618 J = J_OF(IJK)
619 Yi = YG_N(J)
620 ENDIF
621 ENDIF
622
623 IF(DO_K) THEN
624
625
626
627 = X_NODE(4)
628 ya = Y_NODE(4)
629 za = Z_NODE(4)
630
631 N_int_z = 0
632 INTERSECT_Z(IJK) = .FALSE.
633 Q_ID = 1
634 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
635 IF ((INTERSECT_FLAG).AND.(zc/=Zi)) THEN
636 N_int_z = N_int_z + 1
637 INTERSECT_Z(IJK) = .TRUE.
638 zc_backup = Zi
639 Zi = zc
640 ENDIF
641
642 IF(N_int_z /= 1) THEN
643 Zi = UNDEFINED
644 INTERSECT_Z(IJK) = .FALSE.
645 ENDIF
646
647 DO Q_ID = 1, N_POLYGON
648 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Z(IJK),xc,yc,zc)
649 IF(INTERSECT_Z(IJK)) Zi = zc
650 ENDDO
651
652 DO N_USR= 1, N_USR_DEF
653 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Z(IJK),xc,yc,zc)
654 IF(INTERSECT_Z(IJK)) Zi = zc
655 ENDDO
656
657
658
659
660
661
662 IF(TYPE_OF_CELL=='W_MOMENTUM') THEN
663 IF(SNAP(IJK)) THEN
664 INTERSECT_Z(IJK) = .TRUE.
665 K = K_OF(IJK)
666 Zi = ZG_T(K)
667 ENDIF
668 ENDIF
669
670 ENDIF
671
672
673 IF(INTERSECT_X(IJK)) THEN
674 POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
675 POTENTIAL_CUT_CELL_AT(EAST_OF(IJK)) = .TRUE.
676 ENDIF
677
678 IF(INTERSECT_Y(IJK)) THEN
679 POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
680 POTENTIAL_CUT_CELL_AT(NORTH_OF(IJK)) = .TRUE.
681 ENDIF
682
683
684 IF(INTERSECT_Z(IJK)) THEN
685 POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
686 POTENTIAL_CUT_CELL_AT(TOP_OF(IJK)) = .TRUE.
687 ENDIF
688
689
690 RETURN
691
692 END SUBROUTINE INTERSECT
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709 SUBROUTINE CLEAN_INTERSECT(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
710
711 USE param
712 USE param1
713 USE parallel
714 USE constant
715 USE run
716 USE toleranc
717 USE geometry
718 USE indices
719 USE compar
720 USE sendrecv
721 USE quadric
722 USE cutcell
723 USE polygon
724 USE STL
725 USE functions
726
727 IMPLICIT NONE
728 CHARACTER (LEN=*) :: TYPE_OF_CELL
729 INTEGER :: IJK,I,J,K,IM,JM,KM,IP,JP,KP
730 INTEGER :: BCID
731 INTEGER :: IMJK,IPJK,IJMK,IJPK,IJKM,IJKP,IMJPK,IMJKP,IPJMK,IJMKP,IPJKM,IJPKM
732 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb
733 DOUBLE PRECISION :: Xi,Yi,Zi
734 DOUBLE PRECISION :: DFC,DFC_MAX,F4,F6,F7,F8
735 LOGICAL :: CLIP_FLAG,CAD,F_TEST
736
737
738
739
740
741
742 = USE_MSH.OR.USE_STL
743 F_TEST = .TRUE.
744
745
746
747
748
749 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
750
751 I = I_OF(IJK)
752 J = J_OF(IJK)
753 K = K_OF(IJK)
754
755 IM = I - 1
756 JM = J - 1
757 KM = K - 1
758
759 IP = I + 1
760 JP = J + 1
761 KP = K + 1
762
763 IMJK = FUNIJK(IM,J,K)
764 IPJK = FUNIJK(IP,J,K)
765 IJMK = FUNIJK(I,JM,K)
766 IJPK = FUNIJK(I,JP,K)
767 IJKM = FUNIJK(I,J,KM)
768 IJKP = FUNIJK(I,J,KP)
769
770 IMJPK = FUNIJK(IM,JP,K)
771 IMJKP = FUNIJK(IM,J,KP)
772
773 IPJMK = FUNIJK(IP,JM,K)
774 IJMKP = FUNIJK(I,JM,KP)
775
776 IPJKM = FUNIJK(IP,J,KM)
777 IJPKM = FUNIJK(I,JP,KM)
778
779
780 IF(IMJK<1.OR.IMJK>DIMENSION_3) IMJK = IJK
781 IF(IPJK<1.OR.IPJK>DIMENSION_3) IPJK = IJK
782 IF(IJMK<1.OR.IJMK>DIMENSION_3) IJMK = IJK
783 IF(IJPK<1.OR.IJPK>DIMENSION_3) IJPK = IJK
784 IF(IJKM<1.OR.IJKM>DIMENSION_3) IJKM = IJK
785 IF(IJKP<1.OR.IJKP>DIMENSION_3) IJKP = IJK
786
787 IF(IMJPK<1.OR.IMJPK>DIMENSION_3) IMJPK = IJK
788 IF(IMJKP<1.OR.IMJKP>DIMENSION_3) IMJKP = IJK
789
790 IF(IPJMK<1.OR.IPJMK>DIMENSION_3) IPJMK = IJK
791 IF(IJMKP<1.OR.IJMKP>DIMENSION_3) IJMKP = IJK
792
793 IF(IPJKM<1.OR.IPJKM>DIMENSION_3) IPJKM = IJK
794 IF(IJPKM<1.OR.IJPKM>DIMENSION_3) IJPKM = IJK
795
796
797
798
799
800 = X_NODE(7)
801 ya = Y_NODE(7)
802 za = Z_NODE(7)
803
804 xb = X_NODE(8)
805 yb = Y_NODE(8)
806 zb = Z_NODE(8)
807
808 DFC_MAX = TOL_SNAP(1) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
809
810 IF(INTERSECT_X(IJK)) THEN
811
812 DFC = DABS(Xi-xa)
813
814 IF(CAD) F_TEST = (F_AT(IMJK)/=ZERO)
815 IF(DFC < DFC_MAX.AND.F_TEST) THEN
816 IF(PRINT_WARNINGS) THEN
817 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 7'
818 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
819 ENDIF
820
821 INTERSECT_X(IJK) = .FALSE.
822 IF(I>=IMIN1) THEN
823 INTERSECT_X(IMJK) = .FALSE.
824 INTERSECT_Y(IMJK) = .FALSE.
825 INTERSECT_Y(IMJPK) = .FALSE.
826 IF(DO_K) INTERSECT_Z(IMJK) = .FALSE.
827 IF(DO_K) INTERSECT_Z(IMJKP) = .FALSE.
828
829 SNAP(IMJK) = .TRUE.
830 ENDIF
831
832 ENDIF
833
834
835 DFC = DABS(Xi-xb)
836
837 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
838 IF(DFC < DFC_MAX.AND.F_TEST) THEN
839 IF(PRINT_WARNINGS) THEN
840 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 8'
841 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
842 ENDIF
843
844
845 INTERSECT_X(IJK) = .FALSE.
846 INTERSECT_Y(IJK) = .FALSE.
847 IF(DO_K) INTERSECT_Z(IJK) = .FALSE.
848 SNAP(IJK) = .TRUE.
849
850 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
851 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
852 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
853
854
855 ENDIF
856
857 IF(USE_STL.OR.USE_MSH) THEN
858 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,7,F7,CLIP_FLAG,BCID)
859 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
860 IF(F7*F8>TOL_STL**2) INTERSECT_X(IJK) = .FALSE.
861 ENDIF
862
863 ENDIF
864
865
866
867
868
869
870
871
872 = X_NODE(6)
873 ya = Y_NODE(6)
874 za = Z_NODE(6)
875
876 DFC_MAX = TOL_SNAP(2) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
877
878
879 IF(INTERSECT_Y(IJK)) THEN
880
881 DFC = DABS(Yi-ya)
882
883 IF(CAD) F_TEST = (F_AT(IJMK)/=ZERO)
884 IF(DFC < DFC_MAX.AND.F_TEST) THEN
885 IF(PRINT_WARNINGS) THEN
886 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 6'
887 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
888 ENDIF
889
890 INTERSECT_Y(IJK) = .FALSE.
891 IF(J>=JMIN1) THEN
892 INTERSECT_X(IJMK) = .FALSE.
893 INTERSECT_X(IPJMK) = .FALSE.
894 INTERSECT_Y(IJMK) = .FALSE.
895 IF(DO_K) INTERSECT_Z(IJMK) = .FALSE.
896 IF(DO_K) INTERSECT_Z(IJMKP) = .FALSE.
897
898 SNAP(IJMK) = .TRUE.
899 ENDIF
900
901 ENDIF
902
903
904 DFC = DABS(Yi-yb)
905
906 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
907 IF(DFC < DFC_MAX.AND.F_TEST) THEN
908 IF(PRINT_WARNINGS) THEN
909 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 8'
910 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
911 ENDIF
912
913 INTERSECT_X(IJK) = .FALSE.
914 INTERSECT_Y(IJK) = .FALSE.
915 IF(DO_K) INTERSECT_Z(IJK) = .FALSE.
916 SNAP(IJK) = .TRUE.
917
918 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
919 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
920 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
921
922
923 ENDIF
924
925
926 IF(USE_STL.OR.USE_MSH) THEN
927 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,6,F6,CLIP_FLAG,BCID)
928 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
929 IF(F6*F8>TOL_STL**2) INTERSECT_Y(IJK) = .FALSE.
930 ENDIF
931
932 ENDIF
933
934
935 IF(DO_K) THEN
936
937
938
939
940 = X_NODE(4)
941 ya = Y_NODE(4)
942 za = Z_NODE(4)
943
944 DFC_MAX = TOL_SNAP(3) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
945
946 IF(INTERSECT_Z(IJK)) THEN
947
948 DFC = DABS(Zi-Za)
949
950 IF(CAD) F_TEST = (F_AT(IJKM)/=ZERO)
951 IF(DFC < DFC_MAX.AND.F_TEST) THEN
952 IF(PRINT_WARNINGS) THEN
953 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 4'
954 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
955 ENDIF
956
957 INTERSECT_Z(IJK) = .FALSE.
958
959 IF(K>=KMIN1) THEN
960 INTERSECT_X(IJKM) = .FALSE.
961 INTERSECT_X(IPJKM) = .FALSE.
962 INTERSECT_Y(IJKM) = .FALSE.
963 INTERSECT_Y(IJPKM) = .FALSE.
964 INTERSECT_Z(IJKM) = .FALSE.
965
966 SNAP(IJKM) = .TRUE.
967 ENDIF
968
969 ENDIF
970
971
972 DFC = DABS(Zi-Zb)
973
974 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
975 IF(DFC < DFC_MAX.AND.F_TEST) THEN
976 IF(PRINT_WARNINGS) THEN
977 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 8'
978 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
979 ENDIF
980
981 INTERSECT_X(IJK) = .FALSE.
982 INTERSECT_Y(IJK) = .FALSE.
983 INTERSECT_Z(IJK) = .FALSE.
984 SNAP(IJK) = .TRUE.
985
986
987
988
989 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
990 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
991 IF(K<=KMAX1) INTERSECT_Z(IJKP) = .FALSE.
992
993
994 ENDIF
995
996 IF(USE_STL.OR.USE_MSH) THEN
997 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,4,F4,CLIP_FLAG,BCID)
998 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
999 IF(F4*F8>TOL_STL**2) INTERSECT_Z(IJK) = .FALSE.
1000 ENDIF
1001
1002 ENDIF
1003
1004 ENDIF
1005
1006
1007 RETURN
1008
1009 END SUBROUTINE CLEAN_INTERSECT
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026 SUBROUTINE CLEAN_INTERSECT_SCALAR
1027
1028 USE param
1029 USE param1
1030 USE parallel
1031 USE constant
1032 USE run
1033 USE toleranc
1034 USE geometry
1035 USE indices
1036 USE compar
1037 USE sendrecv
1038 USE quadric
1039 USE cutcell
1040 USE polygon
1041 USE STL
1042 USE functions
1043
1044 IMPLICIT NONE
1045 INTEGER :: IJK
1046
1047
1048
1049 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1050 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1051 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1052 call send_recv(Xn_int,2)
1053 call send_recv(Ye_int,2)
1054 call send_recv(Zt_int,2)
1055
1056 DO IJK = IJKSTART3, IJKEND3
1057 CALL SET_SNAP_FLAG(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
1058 END DO
1059
1060 call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1061 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1062 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1063 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1064 call send_recv(Xn_int,2)
1065 call send_recv(Ye_int,2)
1066 call send_recv(Zt_int,2)
1067
1068 DO IJK = IJKSTART3, IJKEND3
1069 CALL REMOVE_INTERSECT_FLAG(IJK)
1070 END DO
1071
1072 call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1073 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1074 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1075 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1076 call send_recv(Xn_int,2)
1077 call send_recv(Ye_int,2)
1078 call send_recv(Zt_int,2)
1079
1080 RETURN
1081
1082 END SUBROUTINE CLEAN_INTERSECT_SCALAR
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098 SUBROUTINE SET_SNAP_FLAG(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
1099
1100 USE param
1101 USE param1
1102 USE parallel
1103 USE constant
1104 USE run
1105 USE toleranc
1106 USE geometry
1107 USE indices
1108 USE compar
1109 USE sendrecv
1110 USE quadric
1111 USE cutcell
1112 USE polygon
1113 USE STL
1114 USE functions
1115
1116 IMPLICIT NONE
1117 CHARACTER (LEN=*) :: TYPE_OF_CELL
1118 INTEGER :: IJK,I,J,K,IM,JM,KM
1119 INTEGER :: BCID
1120 INTEGER :: IMJK,IJMK,IJKM
1121 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb
1122 DOUBLE PRECISION :: Xi,Yi,Zi
1123 DOUBLE PRECISION :: DFC,DFC_MAX,F4,F6,F7,F8
1124 LOGICAL :: CLIP_FLAG,CAD,F_TEST
1125
1126
1127
1128
1129
1130
1131 = USE_MSH.OR.USE_STL
1132 F_TEST = .TRUE.
1133
1134
1135
1136
1137
1138 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1139
1140 I = I_OF(IJK)
1141 J = J_OF(IJK)
1142 K = K_OF(IJK)
1143
1144 IM = I - 1
1145 JM = J - 1
1146 KM = K - 1
1147
1148 IMJK = FUNIJK(IM,J,K)
1149 IJMK = FUNIJK(I,JM,K)
1150 IJKM = FUNIJK(I,J,KM)
1151
1152 IF(IMJK<1.OR.IMJK>DIMENSION_3) IMJK = IJK
1153 IF(IJMK<1.OR.IJMK>DIMENSION_3) IJMK = IJK
1154 IF(IJKM<1.OR.IJKM>DIMENSION_3) IJKM = IJK
1155
1156
1157
1158
1159
1160 = X_NODE(7)
1161 ya = Y_NODE(7)
1162 za = Z_NODE(7)
1163
1164 xb = X_NODE(8)
1165 yb = Y_NODE(8)
1166 zb = Z_NODE(8)
1167
1168
1169 DFC_MAX = TOL_SNAP(1) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1170
1171 IF(INTERSECT_X(IJK)) THEN
1172
1173 DFC = DABS(Xi-xa)
1174
1175 IF(CAD) F_TEST = (F_AT(IMJK)/=ZERO)
1176 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1177 IF(PRINT_WARNINGS) THEN
1178 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 7'
1179 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1180 ENDIF
1181
1182 IF(I>=IMIN1) SNAP(IMJK) = .TRUE.
1183
1184 ENDIF
1185
1186
1187 DFC = DABS(Xi-xb)
1188
1189 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1190 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1191 IF(PRINT_WARNINGS) THEN
1192 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 8'
1193 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1194 ENDIF
1195
1196 SNAP(IJK) = .TRUE.
1197
1198 ENDIF
1199
1200 IF(USE_STL.OR.USE_MSH) THEN
1201 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,7,F7,CLIP_FLAG,BCID)
1202 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1203 IF(F7*F8>TOL_STL**2) INTERSECT_X(IJK) = .FALSE.
1204 ENDIF
1205
1206 ENDIF
1207
1208
1209
1210
1211
1212
1213 = X_NODE(6)
1214 ya = Y_NODE(6)
1215 za = Z_NODE(6)
1216
1217 DFC_MAX = TOL_SNAP(2) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1218
1219
1220 IF(INTERSECT_Y(IJK)) THEN
1221
1222 DFC = DABS(Yi-ya)
1223
1224 IF(CAD) F_TEST = (F_AT(IJMK)/=ZERO)
1225 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1226 IF(PRINT_WARNINGS) THEN
1227 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 6'
1228 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1229 ENDIF
1230
1231 IF(J>=JMIN1) SNAP(IJMK) = .TRUE.
1232
1233 ENDIF
1234
1235
1236 DFC = DABS(Yi-yb)
1237
1238 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1239 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1240 IF(PRINT_WARNINGS) THEN
1241 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 8'
1242 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1243 ENDIF
1244
1245 SNAP(IJK) = .TRUE.
1246
1247
1248 ENDIF
1249
1250
1251 IF(USE_STL.OR.USE_MSH) THEN
1252 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,6,F6,CLIP_FLAG,BCID)
1253 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1254 IF(F6*F8>TOL_STL**2) INTERSECT_Y(IJK) = .FALSE.
1255 ENDIF
1256
1257 ENDIF
1258
1259
1260 IF(DO_K) THEN
1261
1262
1263
1264
1265 = X_NODE(4)
1266 ya = Y_NODE(4)
1267 za = Z_NODE(4)
1268
1269 DFC_MAX = TOL_SNAP(3) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1270
1271 IF(INTERSECT_Z(IJK)) THEN
1272
1273 DFC = DABS(Zi-Za)
1274
1275 IF(CAD) F_TEST = (F_AT(IJKM)/=ZERO)
1276 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1277 IF(PRINT_WARNINGS) THEN
1278 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 4'
1279 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1280 ENDIF
1281
1282 IF(K>=KMIN1) SNAP(IJKM) = .TRUE.
1283
1284 ENDIF
1285
1286 DFC = DABS(Zi-Zb)
1287
1288 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1289 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1290 IF(PRINT_WARNINGS) THEN
1291 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 8'
1292 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1293 ENDIF
1294
1295 SNAP(IJK) = .TRUE.
1296
1297 ENDIF
1298
1299
1300 IF(USE_STL.OR.USE_MSH) THEN
1301 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,4,F4,CLIP_FLAG,BCID)
1302 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1303 IF(F4*F8>TOL_STL**2) INTERSECT_Z(IJK) = .FALSE.
1304 ENDIF
1305
1306 ENDIF
1307
1308 ENDIF
1309
1310
1311 RETURN
1312
1313 END SUBROUTINE SET_SNAP_FLAG
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329 SUBROUTINE REMOVE_INTERSECT_FLAG(IJK)
1330
1331 USE param
1332 USE param1
1333 USE parallel
1334 USE constant
1335 USE run
1336 USE toleranc
1337 USE geometry
1338 USE indices
1339 USE compar
1340 USE sendrecv
1341 USE quadric
1342 USE cutcell
1343 USE polygon
1344 USE STL
1345 USE functions
1346
1347 IMPLICIT NONE
1348 INTEGER :: IJK,I,J,K,IP,JP,KP
1349 INTEGER :: IPJK,IJPK,IJKP
1350
1351
1352 I = I_OF(IJK)
1353 J = J_OF(IJK)
1354 K = K_OF(IJK)
1355
1356 IP = I + 1
1357 JP = J + 1
1358 KP = K + 1
1359
1360 IPJK = FUNIJK(IP,J,K)
1361 IJPK = FUNIJK(I,JP,K)
1362 IJKP = FUNIJK(I,J,KP)
1363
1364 IF(IPJK<1.OR.IPJK>DIMENSION_3) IPJK = IJK
1365 IF(IJPK<1.OR.IJPK>DIMENSION_3) IJPK = IJK
1366 IF(IJKP<1.OR.IJKP>DIMENSION_3) IJKP = IJK
1367
1368 IF(SNAP(IJK)) THEN
1369 INTERSECT_X(IJK) = .FALSE.
1370 INTERSECT_Y(IJK) = .FALSE.
1371 INTERSECT_Z(IJK) = .FALSE.
1372 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
1373 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
1374 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
1375 ENDIF
1376
1377
1378 RETURN
1379
1380 END SUBROUTINE REMOVE_INTERSECT_FLAG
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396 SUBROUTINE OPEN_CUT_CELL_FILES
1397
1398 USE cutcell
1399 USE compar
1400
1401 IMPLICIT NONE
1402
1403 IF(MyPE == PE_IO) THEN
1404 OPEN(UNIT = UNIT_CUT_CELL_LOG, FILE = 'CUT_CELL.LOG')
1405 ENDIF
1406
1407 RETURN
1408
1409
1410 END SUBROUTINE OPEN_CUT_CELL_FILES
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425 SUBROUTINE CLOSE_CUT_CELL_FILES
1426
1427 USE compar
1428 USE cutcell
1429
1430 IMPLICIT NONE
1431
1432 IF(MyPE == PE_IO) THEN
1433 CLOSE(UNIT_CUT_CELL_LOG)
1434 ENDIF
1435
1436 RETURN
1437
1438
1439 END SUBROUTINE CLOSE_CUT_CELL_FILES
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456 SUBROUTINE CAD_INTERSECT(TYPE_OF_CELL,Xint,Yint,Zint)
1457
1458 USE param
1459 USE param1
1460 USE parallel
1461 USE constant
1462 USE run
1463 USE toleranc
1464 USE geometry
1465 USE indices
1466 USE compar
1467 USE sendrecv
1468 USE quadric
1469 USE cutcell
1470 USE polygon
1471 USE stl
1472
1473 USE mpi_utility
1474 USE functions
1475
1476 IMPLICIT NONE
1477 CHARACTER (LEN=*) :: TYPE_OF_CELL
1478 INTEGER :: IJK,I,J,K
1479 INTEGER :: IM,IP,JM,JP,KM,KP,IMJK,IPJK,IJMK,IJPK,IJKM,IJKP
1480 INTEGER :: IJPKP,IPJKP,IPJPK
1481 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb,xc,yc,zc
1482 LOGICAL :: INTERSECT_FLAG,INSIDE_FACET_a,INSIDE_FACET_b
1483
1484 DOUBLE PRECISION :: X1,X2,Y1,Y2,Z1,Z2
1485
1486 DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: Xint,Yint,Zint
1487
1488 INTEGER :: N,I1,I2,J1,J2,K1,K2
1489
1490 DOUBLE PRECISION :: X_OFFSET, Y_OFFSET, Z_OFFSET
1491
1492 DOUBLE PRECISION, DIMENSION(3) :: N4,N6,N7,N8
1493 DOUBLE PRECISION :: CURRENT_F
1494
1495 INTEGER :: N_UNDEFINED, NTOTAL_UNDEFINED,N_PROP
1496 INTEGER, PARAMETER :: N_PROPMAX=1000
1497 LOGICAL:: F_FOUND
1498
1499
1500
1501 = .FALSE.
1502 INTERSECT_Y = .FALSE.
1503 INTERSECT_Z = .FALSE.
1504
1505 Xint = UNDEFINED
1506 Yint = UNDEFINED
1507 Zint = UNDEFINED
1508
1509 F_AT = UNDEFINED
1510
1511 SELECT CASE (TYPE_OF_CELL)
1512 CASE('SCALAR')
1513
1514
1515 X_OFFSET = ZERO
1516 Y_OFFSET = ZERO
1517 Z_OFFSET = ZERO
1518
1519 CASE('U_MOMENTUM')
1520
1521 X_OFFSET = HALF
1522 Y_OFFSET = ZERO
1523 Z_OFFSET = ZERO
1524
1525 CASE('V_MOMENTUM')
1526
1527 X_OFFSET = ZERO
1528 Y_OFFSET = HALF
1529 Z_OFFSET = ZERO
1530
1531 CASE('W_MOMENTUM')
1532
1533 X_OFFSET = ZERO
1534 Y_OFFSET = ZERO
1535 Z_OFFSET = HALF
1536
1537
1538 CASE DEFAULT
1539 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
1540 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1541 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
1542 WRITE(*,*)'SCALAR'
1543 WRITE(*,*)'U_MOMENTUM'
1544 WRITE(*,*)'V_MOMENTUM'
1545 WRITE(*,*)'W_MOMENTUM'
1546 CALL MFIX_EXIT(myPE)
1547 END SELECT
1548
1549
1550 DO N = 1,N_FACETS
1551
1552
1553 X1 = MINVAL(VERTEX(1:3,1,N))
1554 X2 = MAXVAL(VERTEX(1:3,1,N))
1555 Y1 = MINVAL(VERTEX(1:3,2,N))
1556 Y2 = MAXVAL(VERTEX(1:3,2,N))
1557 Z1 = MINVAL(VERTEX(1:3,3,N))
1558 Z2 = MAXVAL(VERTEX(1:3,3,N))
1559
1560
1561 I1 = IEND3
1562 I2 = ISTART3
1563
1564 IF(X2>=ZERO-DX(ISTART3).AND.X1<=XLENGTH+DX(IEND3)) THEN
1565 DO I = ISTART3, IEND3
1566 IP = I+1
1567 IF(XG_E(I)+X_OFFSET*DX(IP)>=X1-TOL_STL) THEN
1568 I1=I
1569 EXIT
1570 ENDIF
1571 ENDDO
1572
1573 DO I = IEND3, ISTART3,-1
1574 IP = I+1
1575 IF(XG_E(I)-DX(I)+X_OFFSET*DX(IP)<=X2+TOL_STL) THEN
1576 I2=I
1577 EXIT
1578 ENDIF
1579 ENDDO
1580 ENDIF
1581
1582
1583 J1 = JEND3
1584 J2 = JSTART3
1585
1586 IF(Y2>=ZERO-DY(JSTART3).AND.Y1<=YLENGTH+DY(JEND3)) THEN
1587 DO J = JSTART3, JEND3
1588 JP = J+1
1589 IF(YG_N(J)+Y_OFFSET*DY(JP)>=Y1-TOL_STL) THEN
1590 J1=J
1591 EXIT
1592 ENDIF
1593 ENDDO
1594
1595 DO J = JEND3, JSTART3,-1
1596 JP=J+1
1597 IF(YG_N(J)-DY(J)+Y_OFFSET*DY(JP)<=Y2+TOL_STL) THEN
1598 J2=J
1599 EXIT
1600 ENDIF
1601 ENDDO
1602 ENDIF
1603
1604 K1 = KEND3
1605 K2 = KSTART3
1606
1607 IF(Z2>=ZERO-DZ(KSTART3).AND.Z1<=ZLENGTH+DZ(KEND3)) THEN
1608 DO K = KSTART3, KEND3
1609 KP=K+1
1610
1611 IF(ZG_T(K)+Z_OFFSET*DZ(KP)>=Z1-TOL_STL) THEN
1612 K1=K
1613 EXIT
1614 ENDIF
1615 ENDDO
1616
1617 DO K = KEND3, KSTART3,-1
1618 KP = K+1
1619 IF(ZG_T(K)-DZ(K)+Z_OFFSET*DZ(KP)<=Z2+TOL_STL) THEN
1620 K2=K
1621 EXIT
1622 ENDIF
1623 ENDDO
1624 ENDIF
1625
1626
1627
1628
1629 DO K=K1,K2
1630 DO J=J1,J2
1631 DO I=I1,I2
1632
1633 IJK = FUNIJK(I,J,K)
1634
1635
1636 IM = MAX0(I - 1 ,ISTART3)
1637 JM = MAX0(J - 1 ,JSTART3)
1638 KM = MAX0(K - 1 ,KSTART3)
1639
1640 IP = MIN0(I + 1 ,IEND3)
1641 JP = MIN0(J + 1 ,JEND3)
1642 KP = MIN0(K + 1 ,KEND3)
1643
1644
1645 IMJK = FUNIJK(IM,J,K)
1646 IPJK = FUNIJK(IP,J,K)
1647 IJMK = FUNIJK(I,JM,K)
1648 IJPK = FUNIJK(I,JP,K)
1649 IJKM = FUNIJK(I,J,KM)
1650 IJKP = FUNIJK(I,J,KP)
1651
1652 IJPKP = FUNIJK(I,JP,KP)
1653 IPJKP = FUNIJK(IP,J,KP)
1654 IPJPK = FUNIJK(IP,JP,K)
1655
1656
1657
1658
1659
1660
1661 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1662
1663
1664
1665
1666 = X_NODE(7)
1667 ya = Y_NODE(7)
1668 za = Z_NODE(7)
1669
1670 xb = X_NODE(8)
1671 yb = Y_NODE(8)
1672 zb = Z_NODE(8)
1673
1674
1675
1676 CALL IS_POINT_INSIDE_FACET(xa,ya,za,N,INSIDE_FACET_a)
1677
1678 IF(INSIDE_FACET_a) THEN
1679
1680 (IMJK) = ZERO
1681
1682 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1683
1684 ENDIF
1685
1686
1687 CALL IS_POINT_INSIDE_FACET(xb,yb,zb,N,INSIDE_FACET_b)
1688
1689 IF(INSIDE_FACET_b) THEN
1690
1691 (IJK) = ZERO
1692
1693 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1694
1695 ENDIF
1696
1697
1698
1699
1700
1701 = .FALSE.
1702
1703 IF(.NOT.(INTERSECT_X(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
1704 CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,N,INTERSECT_FLAG,xc,yc,zc)
1705 ENDIF
1706
1707 IF(INTERSECT_FLAG) THEN
1708 IF(INTERSECT_X(IJK)) THEN
1709 IF(DABS(Xint(IJK)-xc)>TOL_STL) THEN
1710
1711
1712 (IJK) = .FALSE.
1713 ENDIF
1714 ELSE
1715 INTERSECT_X(IJK) = .TRUE.
1716 Xint(IJK) = xc
1717
1718
1719
1720 (1) = xa-xc
1721 N7(2) = ya-yc
1722 N7(3) = za-zc
1723
1724 IF(DABS(F_AT(IMJK))>TOL_F) F_AT(IMJK) = -DOT_PRODUCT(N7,NORM_FACE(:,N))
1725
1726 N8(1) = xb-xc
1727 N8(2) = yb-yc
1728 N8(3) = zb-zc
1729
1730 IF(DABS(F_AT(IJK))>TOL_F) F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,N))
1731
1732
1733 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1734 IF(JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPK,N)
1735 IF(KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJKP,N)
1736 IF(JP<=J2.AND.KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPKP,N)
1737 ENDIF
1738 ENDIF
1739
1740 IF(TYPE_OF_CELL=='U_MOMENTUM') THEN
1741 IF(SNAP(IJK)) THEN
1742 INTERSECT_X(IJK) = .TRUE.
1743 Xn_int(IJK) = XG_E(I)
1744 ENDIF
1745 ENDIF
1746
1747
1748
1749
1750 = X_NODE(6)
1751 ya = Y_NODE(6)
1752 za = Z_NODE(6)
1753
1754
1755
1756 CALL IS_POINT_INSIDE_FACET(xa,ya,za,N,INSIDE_FACET_a)
1757
1758 IF(INSIDE_FACET_a) THEN
1759
1760 (IJMK) = ZERO
1761
1762 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1763
1764 ENDIF
1765
1766
1767
1768
1769
1770
1771
1772 = .FALSE.
1773
1774 IF(.NOT.(INTERSECT_Y(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
1775 CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,N,INTERSECT_FLAG,xc,yc,zc)
1776 ENDIF
1777
1778
1779 IF(INTERSECT_FLAG) THEN
1780
1781 IF(INTERSECT_Y(IJK)) THEN
1782
1783 IF(DABS(Yint(IJK)-yc)>TOL_STL) THEN
1784
1785 INTERSECT_Y(IJK) = .FALSE.
1786
1787 ENDIF
1788
1789 ELSE
1790
1791
1792 INTERSECT_Y(IJK) = .TRUE.
1793 Yint(IJK) = yc
1794
1795
1796
1797 (1) = xa-xc
1798 N6(2) = ya-yc
1799 N6(3) = za-zc
1800
1801 IF(DABS(F_AT(IJMK))>TOL_F) F_AT(IJMK) = -DOT_PRODUCT(N6,NORM_FACE(:,N))
1802
1803 N8(1) = xb-xc
1804 N8(2) = yb-yc
1805 N8(3) = zb-zc
1806
1807 IF(DABS(F_AT(IJK))>TOL_F) F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,N))
1808
1809
1810 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1811 IF(IP<=I2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJK,N)
1812 IF(KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJKP,N)
1813 IF(IP<=I2.AND.KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJKP,N)
1814
1815 ENDIF
1816
1817 ENDIF
1818
1819 IF(TYPE_OF_CELL=='V_MOMENTUM') THEN
1820 IF(SNAP(IJK)) THEN
1821 INTERSECT_Y(IJK) = .TRUE.
1822 Ye_int(IJK) = YG_N(J)
1823 ENDIF
1824 ENDIF
1825
1826 IF(DO_K) THEN
1827
1828
1829
1830 = X_NODE(4)
1831 ya = Y_NODE(4)
1832 za = Z_NODE(4)
1833
1834
1835
1836 CALL IS_POINT_INSIDE_FACET(xa,ya,za,N,INSIDE_FACET_a)
1837
1838 IF(INSIDE_FACET_a) THEN
1839
1840 (IJKM) = ZERO
1841
1842 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1843
1844 ENDIF
1845
1846
1847
1848
1849
1850
1851
1852 = .FALSE.
1853
1854 IF(.NOT.(INTERSECT_Z(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
1855 CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,N,INTERSECT_FLAG,xc,yc,zc)
1856 ENDIF
1857
1858 IF(INTERSECT_FLAG) THEN
1859
1860 IF(INTERSECT_Z(IJK)) THEN
1861
1862 IF(DABS(Zint(IJK)-zc)>TOL_STL) THEN
1863
1864 INTERSECT_Z(IJK) = .FALSE.
1865
1866 ENDIF
1867
1868 ELSE
1869
1870
1871 INTERSECT_Z(IJK) = .TRUE.
1872 Zint(IJK) = zc
1873
1874
1875
1876
1877 (1) = xa-xc
1878 N4(2) = ya-yc
1879 N4(3) = za-zc
1880
1881 IF(DABS(F_AT(IJKM))>TOL_F) F_AT(IJKM) = -DOT_PRODUCT(N4,NORM_FACE(:,N))
1882
1883 N8(1) = xb-xc
1884 N8(2) = yb-yc
1885 N8(3) = zb-zc
1886
1887 IF(DABS(F_AT(IJK))>TOL_F) F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,N))
1888
1889
1890
1891 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,N)
1892 IF(IP<=I2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJK,N)
1893 IF(JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPK,N)
1894 IF(IP<=I2.AND.JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJPK,N)
1895
1896 ENDIF
1897
1898 ENDIF
1899
1900
1901 IF(TYPE_OF_CELL=='W_MOMENTUM') THEN
1902 IF(SNAP(IJK)) THEN
1903 INTERSECT_Z(IJK) = .TRUE.
1904 Zt_int(IJK) = ZG_T(K)
1905 ENDIF
1906 ENDIF
1907
1908 ENDIF
1909
1910
1911 ENDDO
1912 ENDDO
1913 ENDDO
1914
1915
1916 ENDDO
1917
1918 = UNDEFINED
1919
1920
1921
1922 DO IJK = IJKSTART3, IJKEND3
1923 IF(DABS(F_AT(IJK))<TOL_STL) THEN
1924 F_AT(IJK)=ZERO
1925 ENDIF
1926 END DO
1927
1928
1929
1930
1931
1932 call send_recv(F_AT,2)
1933 call send_recv(Xint,2)
1934 call send_recv(Yint,2)
1935 call send_recv(Zint,2)
1936 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1937 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1938 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1939
1940
1941
1942
1943 DO IJK = IJKSTART3, IJKEND3
1944
1945
1946
1947
1948
1949 CALL CLEAN_INTERSECT(IJK,TYPE_OF_CELL,Xint(IJK),Yint(IJK),Zint(IJK))
1950
1951
1952
1953 END DO
1954
1955 call send_recv(F_AT,2)
1956 call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1957 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1958 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1959 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1960
1961 SELECT CASE (CAD_PROPAGATE_ORDER)
1962
1963 CASE (' ')
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981 DO N_PROP=1,N_PROPMAX
1982
1983 DO IJK = IJKSTART3, IJKEND3
1984
1985
1986 = I_OF(IJK)
1987 J = J_OF(IJK)
1988 K = K_OF(IJK)
1989 IF(.NOT.IS_ON_myPE_plus1layer(I,J,K))cycle
1990
1991
1992 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
1993
1994 IMJK = IM_OF(IJK)
1995 IF(F_AT(IMJK)==UNDEFINED.AND.F_AT(IMJK)/=ZERO) F_AT(IMJK)=F_AT(IJK)
1996
1997 IPJK = IP_OF(IJK)
1998 IF(F_AT(IPJK)==UNDEFINED.AND.F_AT(IPJK)/=ZERO) F_AT(IPJK)=F_AT(IJK)
1999
2000 IJMK = JM_OF(IJK)
2001 IF(F_AT(IJMK)==UNDEFINED.AND.F_AT(IJMK)/=ZERO) F_AT(IJMK)=F_AT(IJK)
2002
2003 IJPK = JP_OF(IJK)
2004 IF(F_AT(IJPK)==UNDEFINED.AND.F_AT(IJPK)/=ZERO) F_AT(IJPK)=F_AT(IJK)
2005
2006 IJKM = KM_OF(IJK)
2007 IF(F_AT(IJKM)==UNDEFINED.AND.F_AT(IJKM)/=ZERO) F_AT(IJKM)=F_AT(IJK)
2008
2009 IJKP = KP_OF(IJK)
2010 IF(F_AT(IJKP)==UNDEFINED.AND.F_AT(IJKP)/=ZERO) F_AT(IJKP)=F_AT(IJK)
2011
2012 ENDIF
2013
2014 ENDDO
2015
2016
2017
2018 call send_recv(F_AT,2)
2019
2020
2021
2022 = 0
2023 DO IJK = IJKSTART3, IJKEND3
2024 IF(INTERIOR_CELL_AT(IJK).AND.F_AT(IJK)==UNDEFINED) N_UNDEFINED = N_UNDEFINED + 1
2025 ENDDO
2026
2027 call global_all_sum( N_UNDEFINED, NTOTAL_UNDEFINED )
2028 IF(NTOTAL_UNDEFINED==0) EXIT
2029
2030 ENDDO
2031
2032
2033 call send_recv(F_AT,2)
2034
2035
2036
2037 IF(N_UNDEFINED>0) THEN
2038 WRITE(*,*)'WARNING: UNABLE TO PROPAGATE F_AT ARRAY FROM myPE=.', MyPE
2039 WRITE(*,*)' THIS USUALLY INDICATE A RANK WITH NO ACTIVE CELL'
2040 WRITE(*,*)' YOU MAY NEED TO ADJUST THE GRID PARTITIONNING'
2041 WRITE(*,*)' TO GET BETTER LOAD_BALANCE.'
2042 ENDIF
2043
2044
2045
2046 CASE ('IJK')
2047
2048 DO I=ISTART3,IEND3
2049 DO K=KSTART3,KEND3
2050 F_FOUND = .FALSE.
2051 DO J=JSTART3,JEND3
2052 IJK=FUNIJK(I,J,K)
2053 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2054 F_FOUND = .TRUE.
2055 CURRENT_F = F_AT(IJK)
2056 EXIT
2057 ENDIF
2058 ENDDO
2059 IF(F_FOUND) THEN
2060 DO J=JSTART3,JEND3
2061 IJK=FUNIJK(I,J,K)
2062 IF(F_AT(IJK)==UNDEFINED) THEN
2063 F_AT(IJK)=CURRENT_F
2064 ELSEIF(F_AT(IJK)/=ZERO) THEN
2065 CURRENT_F = F_AT(IJK)
2066 ENDIF
2067 ENDDO
2068 ENDIF
2069 ENDDO
2070 ENDDO
2071
2072 call send_recv(F_AT,2)
2073
2074 DO J=JSTART3,JEND3
2075 DO I=ISTART3,IEND3
2076 F_FOUND = .FALSE.
2077 DO K=KSTART3,KEND3
2078 IJK=FUNIJK(I,J,K)
2079 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2080 F_FOUND = .TRUE.
2081 CURRENT_F = F_AT(IJK)
2082 EXIT
2083 ENDIF
2084 ENDDO
2085 IF(F_FOUND) THEN
2086 DO K=KSTART3,KEND3
2087 IJK=FUNIJK(I,J,K)
2088 IF(F_AT(IJK)==UNDEFINED) THEN
2089 F_AT(IJK)=CURRENT_F
2090 ELSEIF(F_AT(IJK)/=ZERO) THEN
2091 CURRENT_F = F_AT(IJK)
2092 ENDIF
2093 ENDDO
2094 ENDIF
2095 ENDDO
2096 ENDDO
2097
2098 call send_recv(F_AT,2)
2099
2100 DO K=KSTART3,KEND3
2101 DO J=JSTART3,JEND3
2102 F_FOUND = .FALSE.
2103 DO I=ISTART3,IEND3
2104 IJK=FUNIJK(I,J,K)
2105 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2106 F_FOUND = .TRUE.
2107 CURRENT_F = F_AT(IJK)
2108 EXIT
2109 ENDIF
2110 ENDDO
2111 IF(F_FOUND) THEN
2112 DO I=ISTART3,IEND3
2113 IJK=FUNIJK(I,J,K)
2114 IF(F_AT(IJK)==UNDEFINED) THEN
2115 F_AT(IJK)=CURRENT_F
2116 ELSEIF(F_AT(IJK)/=ZERO) THEN
2117 CURRENT_F = F_AT(IJK)
2118 ENDIF
2119 ENDDO
2120 ENDIF
2121 ENDDO
2122 ENDDO
2123
2124 call send_recv(F_AT,2)
2125
2126 CASE ('JKI')
2127
2128 DO J=JSTART3,JEND3
2129 DO I=ISTART3,IEND3
2130 F_FOUND = .FALSE.
2131 DO K=KSTART3,KEND3
2132 IJK=FUNIJK(I,J,K)
2133 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2134 F_FOUND = .TRUE.
2135 CURRENT_F = F_AT(IJK)
2136 EXIT
2137 ENDIF
2138 ENDDO
2139 IF(F_FOUND) THEN
2140 DO K=KSTART3,KEND3
2141 IJK=FUNIJK(I,J,K)
2142 IF(F_AT(IJK)==UNDEFINED) THEN
2143 F_AT(IJK)=CURRENT_F
2144 ELSEIF(F_AT(IJK)/=ZERO) THEN
2145 CURRENT_F = F_AT(IJK)
2146 ENDIF
2147 ENDDO
2148 ENDIF
2149 ENDDO
2150 ENDDO
2151
2152 call send_recv(F_AT,2)
2153
2154 DO K=KSTART3,KEND3
2155 DO J=JSTART3,JEND3
2156 F_FOUND = .FALSE.
2157 DO I=ISTART3,IEND3
2158 IJK=FUNIJK(I,J,K)
2159 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2160 F_FOUND = .TRUE.
2161 CURRENT_F = F_AT(IJK)
2162 EXIT
2163 ENDIF
2164 ENDDO
2165 IF(F_FOUND) THEN
2166 DO I=ISTART3,IEND3
2167 IJK=FUNIJK(I,J,K)
2168 IF(F_AT(IJK)==UNDEFINED) THEN
2169 F_AT(IJK)=CURRENT_F
2170 ELSEIF(F_AT(IJK)/=ZERO) THEN
2171 CURRENT_F = F_AT(IJK)
2172 ENDIF
2173 ENDDO
2174 ENDIF
2175 ENDDO
2176 ENDDO
2177
2178 call send_recv(F_AT,2)
2179
2180 DO I=ISTART3,IEND3
2181 DO K=KSTART3,KEND3
2182 F_FOUND = .FALSE.
2183 DO J=JSTART3,JEND3
2184 IJK=FUNIJK(I,J,K)
2185 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2186 F_FOUND = .TRUE.
2187 CURRENT_F = F_AT(IJK)
2188 EXIT
2189 ENDIF
2190 ENDDO
2191 IF(F_FOUND) THEN
2192 DO J=JSTART3,JEND3
2193 IJK=FUNIJK(I,J,K)
2194 IF(F_AT(IJK)==UNDEFINED) THEN
2195 F_AT(IJK)=CURRENT_F
2196 ELSEIF(F_AT(IJK)/=ZERO) THEN
2197 CURRENT_F = F_AT(IJK)
2198 ENDIF
2199 ENDDO
2200 ENDIF
2201 ENDDO
2202 ENDDO
2203
2204 call send_recv(F_AT,2)
2205
2206
2207 CASE ('KIJ')
2208
2209
2210
2211 DO K=KSTART3,KEND3
2212 DO J=JSTART3,JEND3
2213 F_FOUND = .FALSE.
2214 DO I=ISTART3,IEND3
2215 IJK=FUNIJK(I,J,K)
2216 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2217 F_FOUND = .TRUE.
2218 CURRENT_F = F_AT(IJK)
2219 EXIT
2220 ENDIF
2221 ENDDO
2222 IF(F_FOUND) THEN
2223 DO I=ISTART3,IEND3
2224 IJK=FUNIJK(I,J,K)
2225 IF(F_AT(IJK)==UNDEFINED) THEN
2226 F_AT(IJK)=CURRENT_F
2227 ELSEIF(F_AT(IJK)/=ZERO) THEN
2228 CURRENT_F = F_AT(IJK)
2229 ENDIF
2230 ENDDO
2231 ENDIF
2232 ENDDO
2233 ENDDO
2234
2235 call send_recv(F_AT,2)
2236
2237 DO I=ISTART3,IEND3
2238 DO K=KSTART3,KEND3
2239 F_FOUND = .FALSE.
2240 DO J=JSTART3,JEND3
2241 IJK=FUNIJK(I,J,K)
2242 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2243 F_FOUND = .TRUE.
2244 CURRENT_F = F_AT(IJK)
2245 EXIT
2246 ENDIF
2247 ENDDO
2248 IF(F_FOUND) THEN
2249 DO J=JSTART3,JEND3
2250 IJK=FUNIJK(I,J,K)
2251 IF(F_AT(IJK)==UNDEFINED) THEN
2252 F_AT(IJK)=CURRENT_F
2253 ELSEIF(F_AT(IJK)/=ZERO) THEN
2254 CURRENT_F = F_AT(IJK)
2255 ENDIF
2256 ENDDO
2257 ENDIF
2258 ENDDO
2259 ENDDO
2260
2261
2262 call send_recv(F_AT,2)
2263
2264 DO J=JSTART3,JEND3
2265 DO I=ISTART3,IEND3
2266 F_FOUND = .FALSE.
2267 DO K=KSTART3,KEND3
2268 IJK=FUNIJK(I,J,K)
2269 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2270 F_FOUND = .TRUE.
2271 CURRENT_F = F_AT(IJK)
2272 EXIT
2273 ENDIF
2274 ENDDO
2275 IF(F_FOUND) THEN
2276 DO K=KSTART3,KEND3
2277 IJK=FUNIJK(I,J,K)
2278 IF(F_AT(IJK)==UNDEFINED) THEN
2279 F_AT(IJK)=CURRENT_F
2280 ELSEIF(F_AT(IJK)/=ZERO) THEN
2281 CURRENT_F = F_AT(IJK)
2282 ENDIF
2283 ENDDO
2284 ENDIF
2285 ENDDO
2286 ENDDO
2287
2288
2289 call send_recv(F_AT,2)
2290
2291 CASE DEFAULT
2292 IF(myPE == PE_IO) THEN
2293 WRITE(*,*)'CAD_INTERSECT.'
2294 WRITE(*,*)'UNKNOWN CAD_PROPAGATE_ORDER:',CAD_PROPAGATE_ORDER
2295 WRITE(*,*)'ACCEPTABLE VALUES:'
2296 WRITE(*,*)'IJK'
2297 WRITE(*,*)'JKI'
2298 WRITE(*,*)'KIJ'
2299 ENDIF
2300
2301
2302 END SELECT
2303
2304
2305
2306 call send_recv(F_AT,2)
2307
2308 RETURN
2309
2310 END SUBROUTINE CAD_INTERSECT
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326 SUBROUTINE ADD_FACET_AND_SET_BC_ID(IJK,N)
2327
2328 USE param
2329 USE param1
2330 USE parallel
2331 USE constant
2332 USE run
2333 USE toleranc
2334 USE geometry
2335 USE indices
2336 USE compar
2337 USE sendrecv
2338 USE quadric
2339 USE cutcell
2340 USE polygon
2341 USE stl
2342 USE mpi_utility
2343
2344 IMPLICIT NONE
2345 INTEGER :: IJK,N
2346
2347
2348 BC_ID(IJK) = BC_ID_STL_FACE(N)
2349
2350 IF(N_FACET_AT(IJK)<DIM_FACETS_PER_CELL) THEN
2351
2352 N_FACET_AT(IJK) = N_FACET_AT(IJK) + 1
2353 LIST_FACET_AT(IJK,N_FACET_AT(IJK)) = N
2354
2355 ELSE
2356
2357 WRITE(*,*) ' FATAL ERROR: TOO MANY FACETS IN CELL: ', IJK
2358 CALL MFIX_EXIT(myPE)
2359
2360 ENDIF
2361
2362
2363
2364 END SUBROUTINE ADD_FACET_AND_SET_BC_ID
2365
2366
2367