File: N:\mfix\model\cartesian_grid\cut_cell_preprocessing.f
1 MODULE CUT_CELL_PREPROC
2 USE exit, ONLY: mfix_exit
3 CONTAINS
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21 SUBROUTINE CUT_CELL_PREPROCESSING
22
23 USE cdist
24 USE compar
25 USE constant
26 USE cutcell
27 USE fldvar
28 USE functions
29 USE geometry
30 USE indices
31 USE parallel
32 USE param
33 USE param1
34 USE quadric
35 USE run
36 USE sendrecv
37 USE toleranc
38 USE vtk
39
40 IMPLICIT NONE
41
42 INTEGER :: SAFE_MODE_COUNT
43 DOUBLE PRECISION :: CPU_PP_START,CPU_PP_END
44
45 IF(.NOT.CG_HEADER_WAS_PRINTED) CALL PRINT_CG_HEADER
46
47 CALL CPU_TIME (CPU_PP_START)
48
49 CALL OPEN_CUT_CELL_FILES
50
51 CALL ALLOCATE_CUT_CELL_ARRAYS
52
53 CALL DEFINE_QUADRICS
54
55 CALL SET_3D_CUT_CELL_FLAGS
56
57 CALL GATHER_DATA
58
59
60
61
62
63 IF(WRITE_VTK_FILES.AND.(.NOT.BDIST_IO)) THEN
64 CALL WRITE_CUT_SURFACE_VTK
65 ENDIF
66
67 CALL SET_3D_CUT_U_CELL_FLAGS
68 CALL SET_3D_CUT_V_CELL_FLAGS
69 IF(DO_K) CALL SET_3D_CUT_W_CELL_FLAGS
70
71 CALL SET_3D_CUT_CELL_TREATMENT_FLAGS
72
73 CALL GET_3D_ALPHA_U_CUT_CELL
74 CALL GET_3D_ALPHA_V_CUT_CELL
75 IF(DO_K) CALL GET_3D_ALPHA_W_CUT_CELL
76
77 CALL SET_GHOST_CELL_FLAGS
78
79 CALL SET_ODXYZ_U_CUT_CELL
80 CALL SET_ODXYZ_V_CUT_CELL
81 IF(DO_K) CALL SET_ODXYZ_W_CUT_CELL
82
83 CALL GET_U_MASTER_CELLS
84 CALL GET_V_MASTER_CELLS
85 IF(DO_K) CALL GET_W_MASTER_CELLS
86
87 CALL SEND_RECEIVE_CUT_CELL_VARIABLES
88
89 CALL GET_DISTANCE_TO_WALL
90
91 CALL PRINT_GRID_STATISTICS
92
93 CALL CG_GET_BC_AREA
94
95 CALL FLOW_TO_VEL(.FALSE.)
96
97 CALL CG_FLOW_TO_VEL
98
99 CALL CONVERT_CG_MI_TO_PS
100
101 CALL CPU_TIME (CPU_PP_END)
102
103 IF(myPE == PE_IO) THEN
104 WRITE(*,20)'CARTESIAN GRID PRE-PROCESSING COMPLETED IN ',CPU_PP_END - CPU_PP_START, ' SECONDS.'
105 WRITE(*,10)'============================================================================='
106 ENDIF
107
108 IF(myPE == PE_IO) THEN
109
110 SAFE_MODE_COUNT = SUM(CG_SAFE_MODE)
111
112 IF(SAFE_MODE_COUNT>0) THEN
113
114
115 WRITE(*,10)'######################################################################'
116 WRITE(*,10)'######################################################################'
117 WRITE(*,10)'## ##'
118 WRITE(*,10)'## || ##'
119 WRITE(*,10)'## || ##'
120 WRITE(*,10)'## \/ ##'
121 WRITE(*,10)'## ##'
122 WRITE(*,10)'## ===> WARNING: RUNNING CARTESIAN GRID IN SAFE MODE ! <=== ##'
123 WRITE(*,10)'## ##'
124 WRITE(*,10)'## SAFE MODE ACTIVATED FOR : ##'
125 IF(CG_SAFE_MODE(1)==1) WRITE(*,10)'## - All scalar quantities ##'
126 IF(CG_SAFE_MODE(3)==1) WRITE(*,10)'## - X-Velocity (Gas and Solids) ##'
127 IF(CG_SAFE_MODE(4)==1) WRITE(*,10)'## - Y-Velocity (Gas and Solids) ##'
128 IF(CG_SAFE_MODE(5)==1) WRITE(*,10)'## - Z-Velocity (Gas and Solids) ##'
129 WRITE(*,10)'## ##'
130 WRITE(*,10)'## /\ ##'
131 WRITE(*,10)'## || ##'
132 WRITE(*,10)'## || ##'
133 WRITE(*,10)'## ##'
134 WRITE(*,10)'######################################################################'
135 WRITE(*,10)'######################################################################'
136
137 ENDIF
138 ENDIF
139
140 RETURN
141
142 10 FORMAT(1X,A)
143 20 FORMAT(1X,A,F8.2,A)
144
145 END SUBROUTINE CUT_CELL_PREPROCESSING
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161 SUBROUTINE CG_FLOW_TO_VEL
162
163 USE bc
164 USE compar
165 USE constant
166 USE cutcell
167 USE eos, ONLY: EOSG
168 USE fldvar
169 USE functions
170 USE funits
171 USE geometry
172 USE indices
173 USE mpi_utility
174 USE parallel
175 USE param
176 USE param1
177 USE physprop
178 USE quadric
179 USE run
180 USE scales
181 USE sendrecv
182 USE toleranc
183 USE vtk
184
185 IMPLICIT NONE
186
187
188
189
190
191
192
193
194
195
196
197 INTEGER :: M, BCV
198
199 DOUBLE PRECISION :: VOLFLOW
200
201 DOUBLE PRECISION :: EPS
202
203 DOUBLE PRECISION :: MW
204
205
206
207 DO BCV = 1, DIMENSION_BC
208
209 IF (BC_TYPE_ENUM(BCV)==CG_MI) THEN
210
211 IF(BC_VELMAG_g(BCV)==UNDEFINED) THEN
212
213
214
215 IF (BC_MASSFLOW_G(BCV) /= UNDEFINED) THEN
216 IF (RO_G0 /= UNDEFINED) THEN
217 VOLFLOW = BC_MASSFLOW_G(BCV)/RO_G0
218 ELSE
219 IF (BC_P_G(BCV)/=UNDEFINED .AND. BC_T_G(BCV)/=UNDEFINED) &
220 THEN
221 IF (MW_AVG == UNDEFINED) THEN
222 MW = CALC_MW(BC_X_G,DIMENSION_BC,BCV,NMAX(0),MW_G)
223 ELSE
224 MW = MW_AVG
225 ENDIF
226 VOLFLOW = BC_MASSFLOW_G(BCV)/EOSG(MW,(BC_P_G(BCV)-P_REF), &
227 BC_T_G(BCV))
228 ELSE
229 IF (BC_TYPE_ENUM(BCV) == CG_MO) THEN
230 IF (BC_MASSFLOW_G(BCV) == ZERO) THEN
231 VOLFLOW = ZERO
232 ENDIF
233 ELSE
234 IF(DMP_LOG)WRITE (UNIT_LOG, 1020) BCV
235 call mfix_exit(myPE)
236 ENDIF
237 ENDIF
238 ENDIF
239
240
241
242 IF (BC_VOLFLOW_G(BCV) /= UNDEFINED) THEN
243 IF (.NOT.COMPARE(VOLFLOW,BC_VOLFLOW_G(BCV))) THEN
244 IF(DMP_LOG)WRITE (UNIT_LOG, 1000) BCV, VOLFLOW, BC_VOLFLOW_G(BCV)
245 call mfix_exit(myPE)
246 ENDIF
247 ELSE
248 BC_VOLFLOW_G(BCV) = VOLFLOW
249 ENDIF
250 ENDIF
251
252
253
254 IF (BC_VOLFLOW_G(BCV) /= UNDEFINED) THEN
255 IF (BC_EP_G(BCV) /= UNDEFINED) THEN
256 BC_VELMAG_g(BCV) = BC_VOLFLOW_G(BCV)/(BC_AREA(BCV)*BC_EP_G(BCV))
257 ELSE
258 RETURN
259 ENDIF
260 ENDIF
261
262 ENDIF
263
264
265
266
267
268 DO M = 1, MMAX
269
270 IF(BC_VELMAG_s(BCV,M)==UNDEFINED) THEN
271
272
273
274 IF (BC_MASSFLOW_S(BCV,M) /= UNDEFINED) THEN
275 IF (RO_S0(M) /= UNDEFINED) THEN
276 VOLFLOW = BC_MASSFLOW_S(BCV,M)/RO_S0(M)
277 ELSE
278 RETURN
279 ENDIF
280
281
282
283 IF (BC_VOLFLOW_S(BCV,M) /= UNDEFINED) THEN
284 IF (.NOT.COMPARE(VOLFLOW,BC_VOLFLOW_S(BCV,M))) THEN
285 IF(DMP_LOG)WRITE(UNIT_LOG,1200)BCV,VOLFLOW,M,BC_VOLFLOW_S(BCV,M)
286 call mfix_exit(myPE)
287 ENDIF
288 ELSE
289 BC_VOLFLOW_S(BCV,M) = VOLFLOW
290 ENDIF
291 ENDIF
292
293 IF (BC_ROP_S(BCV,M)==UNDEFINED .AND. MMAX==1) BC_ROP_S(BCV,M)&
294 = (ONE - BC_EP_G(BCV))*RO_S0(M)
295 IF (BC_VOLFLOW_S(BCV,M) /= UNDEFINED) THEN
296 IF (BC_ROP_S(BCV,M) /= UNDEFINED) THEN
297 EPS = BC_ROP_S(BCV,M)/RO_S0(M)
298 IF (EPS /= ZERO) THEN
299 BC_VELMAG_s(BCV,M) = BC_VOLFLOW_S(BCV,M)/(BC_AREA(BCV)*EPS)
300 ELSE
301 IF (BC_VOLFLOW_S(BCV,M) == ZERO) THEN
302 BC_VELMAG_s(BCV,M) = ZERO
303 ELSE
304 IF(DMP_LOG)WRITE (UNIT_LOG, 1250) BCV, M
305 call mfix_exit(myPE)
306 ENDIF
307 ENDIF
308 ELSE
309 IF (BC_VOLFLOW_S(BCV,M) == ZERO) THEN
310 BC_VELMAG_s(BCV,M) = ZERO
311 ELSE
312 IF(DMP_LOG)WRITE (UNIT_LOG, 1260) BCV, M
313 call mfix_exit(myPE)
314 ENDIF
315 ENDIF
316 ENDIF
317
318 ENDIF
319 END DO
320 ENDIF
321 END DO
322
323 100 FORMAT(1X,A,I8)
324 110 FORMAT(1X,A,A)
325 120 FORMAT(1X,A,F14.8,/)
326 130 FORMAT(1X,A,I8,F14.8,/)
327
328 1000 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
329 ' Computed volumetric flow is not equal to specified value',/,&
330 ' Value computed from mass flow = ',G14.7,/,&
331 ' Specified value (BC_VOLFLOW_g) = ',G14.7,/1X,70('*')/)
332
333
334 1020 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,&
335 ' BC_P_g, BC_T_g, and BC_X_g',/' should be specified',/1X,70('*')/)
336
337
338 1200 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
339 ' Computed volumetric flow is not equal to specified value',/,&
340 ' Value computed from mass flow = ',G14.7,/,&
341 ' Specified value (BC_VOLFLOW_s',I1,') = ',G14.7,/1X,70('*')/)
342
343 1250 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
344 ' Non-zero vol. or mass flow specified with BC_ROP_s',&
345 I1,' = 0.',/1X,70('*')/)
346 1260 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
347 ' BC_ROP_s',I1,' not specified',/1X,70('*')/)
348 RETURN
349
350 END SUBROUTINE CG_FLOW_TO_VEL
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366 SUBROUTINE CONVERT_CG_MI_TO_PS
367
368 USE bc
369 USE compar
370 USE constant
371 USE cutcell
372 USE eos, only: EOSG
373 USE fldvar
374 USE functions
375 USE funits
376 USE geometry
377 USE indices
378 USE mpi_utility
379 USE parallel
380 USE param
381 USE param1
382 USE physprop
383 USE ps
384 USE quadric
385 USE run
386 USE scales
387 USE sendrecv
388 USE toleranc
389 USE vtk
390
391 IMPLICIT NONE
392
393
394
395
396
397
398
399
400
401
402
403 INTEGER :: IJK, M, NN, BCV
404
405 INTEGER :: iproc
406 INTEGER :: NPS,PSV
407
408
409
410
411
412 #ifdef MPI
413 CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
414 #endif
415
416
417
418 do iproc = 0,NumPEs-1
419 if (MyPE==iproc) Then
420
421
422
423 = 0
424
425 PS_LP: do PSV = 1, DIMENSION_PS
426 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
427 NPS = PSV
428 enddo PS_LP
429
430
431
432
433
434 DO IJK = ijkstart3, ijkend3
435 BCV = BC_ID(IJK)
436 IF(BCV>0) THEN
437 IF(CG_MI_CONVERTED_TO_PS(BCV).AND.INTERIOR_CELL_AT(IJK).AND.VOL(IJK)>ZERO) THEN
438
439 NPS = NPS + 1
440
441
442
443 (NPS) = .TRUE.
444
445 POINT_SOURCE = .TRUE.
446
447 PS_I_w(NPS) = I_OF(IJK)
448 PS_I_e(NPS) = I_OF(IJK)
449 PS_J_s(NPS) = J_OF(IJK)
450 PS_J_n(NPS) = J_OF(IJK)
451 PS_K_b(NPS) = K_OF(IJK)
452 PS_K_t(NPS) = K_OF(IJK)
453
454 PS_VOLUME(NPS) = VOL(IJK)
455
456 PS_MASSFLOW_g(NPS) = BC_MASSFLOW_g(BCV) * VOL(IJK) / BC_VOL(BCV)
457
458 PS_T_g(NPS) = BC_T_g(BCV)
459
460 IF(BC_U_g(BCV)==UNDEFINED) THEN
461 PS_U_g(NPS) = Normal_S(IJK,1)
462 ELSE
463 PS_U_g(NPS) = BC_U_g(BCV)
464 ENDIF
465
466 IF(BC_V_g(BCV)==UNDEFINED) THEN
467 PS_V_g(NPS) = Normal_S(IJK,2)
468 ELSE
469 PS_V_g(NPS) = BC_V_g(BCV)
470 ENDIF
471
472 IF(BC_W_g(BCV)==UNDEFINED) THEN
473 PS_W_g(NPS) = Normal_S(IJK,3)
474 ELSE
475 PS_W_g(NPS) = BC_W_g(BCV)
476 ENDIF
477
478 DO NN=1,NMAX(0)
479 PS_X_g(NPS,NN) = BC_X_g(BCV,NN)
480 ENDDO
481
482 DO M=1, MMAX
483 PS_MASSFLOW_s(NPS,M) = BC_MASSFLOW_s(BCV,M) * VOL(IJK) / BC_VOL(BCV)
484
485 PS_T_s(NPS,1) = BC_T_s(BCV,M)
486
487 IF(BC_U_s(BCV,M)==UNDEFINED) THEN
488 PS_U_s(NPS,M) = Normal_S(IJK,1)
489 ELSE
490 PS_U_s(NPS,M) = BC_U_s(BCV,M)
491 ENDIF
492
493 IF(BC_V_s(BCV,M)==UNDEFINED) THEN
494 PS_V_s(NPS,M) = Normal_S(IJK,2)
495 ELSE
496 PS_V_s(NPS,M) = BC_V_s(BCV,M)
497 ENDIF
498
499 IF(BC_W_s(BCV,M)==UNDEFINED) THEN
500 PS_W_s(NPS,M) = Normal_S(IJK,3)
501 ELSE
502 PS_W_s(NPS,M) = BC_W_s(BCV,M)
503 ENDIF
504
505
506 DO NN=1,NMAX(M)
507 PS_X_s(NPS,M,NN) = BC_X_s(BCV,M,NN)
508 ENDDO
509
510 ENDDO
511
512
513 ENDIF
514 ENDIF
515
516 ENDDO
517
518 endif
519
520 #ifdef MPI
521 CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
522 #endif
523 call bcast(POINT_SOURCE,iproc)
524 call bcast(PS_DEFINED,iproc)
525 call bcast(PS_I_w,iproc)
526 call bcast(PS_I_e,iproc)
527 call bcast(PS_J_s,iproc)
528 call bcast(PS_J_n,iproc)
529 call bcast(PS_K_b,iproc)
530 call bcast(PS_K_t,iproc)
531 call bcast(PS_MASSFLOW_g,iproc)
532 call bcast(PS_U_g,iproc)
533 call bcast(PS_V_g,iproc)
534 call bcast(PS_W_g,iproc)
535 call bcast(PS_X_g,iproc)
536 call bcast(PS_T_g,iproc)
537 call bcast(PS_MASSFLOW_s,iproc)
538 call bcast(PS_U_s,iproc)
539 call bcast(PS_V_s,iproc)
540 call bcast(PS_W_s,iproc)
541 call bcast(PS_X_s,iproc)
542 call bcast(PS_T_s,iproc)
543 call bcast(PS_VOLUME,iproc)
544
545 enddo
546
547 #ifdef MPI
548 CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
549 #endif
550
551
552
553 RETURN
554
555 END SUBROUTINE CONVERT_CG_MI_TO_PS
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571 SUBROUTINE CONVERT_CG_MI_TO_PS_PE
572
573 USE bc
574 USE compar
575 USE constant
576 USE cutcell
577 USE eos, ONLY: EOSG
578 USE fldvar
579 USE functions
580 USE funits
581 USE geometry
582 USE indices
583 USE mpi_utility
584 USE parallel
585 USE param
586 USE param1
587 USE physprop
588 USE ps
589 USE quadric
590 USE run
591 USE scales
592 USE sendrecv
593 USE toleranc
594 USE vtk
595
596 IMPLICIT NONE
597
598
599
600
601
602
603
604
605
606
607
608 INTEGER :: IJK, BCV
609
610 INTEGER :: NPS,PSV
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627 = 0
628
629 PS_LP: do PSV = 1, DIMENSION_PS
630 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
631 NPS = PSV
632 enddo PS_LP
633
634 print *,'Last PS=',NPS
635
636
637
638
639
640 DO IJK = ijkstart3, ijkend3
641 BCV = BC_ID(IJK)
642 IF(BCV>0) THEN
643 IF(CG_MI_CONVERTED_TO_PS(BCV).AND.INTERIOR_CELL_AT(IJK).AND.VOL(IJK)>ZERO) THEN
644
645 NPS = NPS + 1
646
647 PS_DEFINED(NPS) = .TRUE.
648
649 POINT_SOURCE = .TRUE.
650
651 PS_I_w(NPS) = I_OF(IJK)
652 PS_I_e(NPS) = I_OF(IJK)
653 PS_J_s(NPS) = J_OF(IJK)
654 PS_J_n(NPS) = J_OF(IJK)
655 PS_K_b(NPS) = K_OF(IJK)
656 PS_K_t(NPS) = K_OF(IJK)
657
658 PS_MASSFLOW_g(NPS) = BC_MASSFLOW_g(BCV) * VOL(IJK) / BC_VOL(BCV)
659
660 PS_VOLUME(NPS) = VOL(IJK)
661
662 PS_T_g(NPS) = BC_T_g(BCV)
663
664 IF(BC_U_g(NPS)==UNDEFINED) THEN
665 PS_U_g(NPS) = Normal_S(IJK,1)
666 ELSE
667 PS_U_g(NPS) = BC_U_g(NPS)
668 ENDIF
669
670 IF(BC_V_g(NPS)==UNDEFINED) THEN
671 PS_V_g(NPS) = Normal_S(IJK,2)
672 ELSE
673 PS_V_g(NPS) = BC_V_g(NPS)
674 ENDIF
675
676 IF(BC_W_g(NPS)==UNDEFINED) THEN
677 PS_W_g(NPS) = Normal_S(IJK,3)
678 ELSE
679 PS_W_g(NPS) = BC_W_g(NPS)
680 ENDIF
681
682
683 (NPS,1) = 0.0
684
685 PS_T_s(NPS,1) = 298.0
686
687 PS_U_s(NPS,1) = Normal_S(IJK,1)
688 PS_V_s(NPS,1) = Normal_S(IJK,2)
689 PS_W_s(NPS,1) = Normal_S(IJK,3)
690
691 PS_U_s(NPS,1) = ZERO
692 PS_V_s(NPS,1) = ZERO
693 PS_W_s(NPS,1) = ZERO
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720 print*,'PS created:',NPS,PS_MASSFLOW_g(NPS),PS_VOLUME(NPS),PS_I_w(NPS),PS_J_n(NPS),PS_K_b(NPS), &
721 INTERIOR_CELL_AT(IJK),PS_U_g(NPS),PS_V_g(NPS),PS_W_g(NPS), &
722 CUT_CELL_AT(IJK),Normal_S(IJK,1),Normal_S(IJK,2),Normal_S(IJK,3)
723
724
725
726 ENDIF
727 ENDIF
728 ENDDO
729
730
731
732
733
734
735
736
737 FORMAT(1X,A,I8)
738 110 FORMAT(1X,A,A)
739 120 FORMAT(1X,A,F14.8,/)
740 130 FORMAT(1X,A,I8,F14.8,/)
741
742
743 1000 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
744 ' Computed volumetric flow is not equal to specified value',/,&
745 ' Value computed from mass flow = ',G14.7,/,&
746 ' Specified value (BC_VOLFLOW_g) = ',G14.7,/1X,70('*')/)
747
748
749 1020 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,&
750 ' BC_P_g, BC_T_g, and BC_X_g',/' should be specified',/1X,70('*')/)
751
752
753 1200 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
754 ' Computed volumetric flow is not equal to specified value',/,&
755 ' Value computed from mass flow = ',G14.7,/,&
756 ' Specified value (BC_VOLFLOW_s',I1,') = ',G14.7,/1X,70('*')/)
757
758 1250 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
759 ' Non-zero vol. or mass flow specified with BC_ROP_s',&
760 I1,' = 0.',/1X,70('*')/)
761 1260 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
762 ' BC_ROP_s',I1,' not specified',/1X,70('*')/)
763 RETURN
764
765 END SUBROUTINE CONVERT_CG_MI_TO_PS_PE
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780 SUBROUTINE PRINT_CG_HEADER
781
782 USE param
783 USE param1
784 USE parallel
785 USE compar
786 USE cutcell
787
788 IMPLICIT NONE
789
790
791
792
793 IF(myPE == PE_IO) THEN
794
795 WRITE(*,10)'============================================================================='
796 WRITE(*,10)' ____ ___________ _ ____ ____ __________ __________ '
797 WRITE(*,10)' | \ / | (_) \ \ / / | | | | '
798 WRITE(*,10)' | \ / ______| ___ \ \ / / | ______| | ______| '
799 WRITE(*,10)' | \ / |______ | | \ \/ / | | | | '
800 WRITE(*,10)' | \ / | | | \ / === | | | | ____ '
801 WRITE(*,10)' | |\ \/ /| ______| | | / \ | | | | |_ | '
802 WRITE(*,10)' | | \ / | | | | / /\ \ | |______ | |___| | '
803 WRITE(*,10)' | | \ / | | | | / / \ \ | | | | '
804 WRITE(*,10)' |___| \__/ |___| |___| /___/ \___\ |__________| |__________| '
805 WRITE(*,10)' '
806 WRITE(*,10)'============================================================================='
807 WRITE(*,10)'MFIX WITH CARTESIAN GRID IMPLEMENTATION.'
808
809
810 IF(RE_INDEXING) THEN
811 WRITE(*,10)'RE-INDEXING IS TURNED ON.'
812
813
814
815
816
817
818 ELSE
819 WRITE(*,10)'RE-INDEXING IS TURNED OFF.'
820 ENDIF
821
822 ENDIF
823
824 CG_HEADER_WAS_PRINTED = .TRUE.
825
826 10 FORMAT(1X,A)
827
828 RETURN
829
830 END SUBROUTINE PRINT_CG_HEADER
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854 SUBROUTINE EVAL_F(METHOD,x1,x2,x3,Q,f,CLIP_FLAG)
855
856 USE compar
857 USE cutcell
858 USE parallel
859 USE param1, only: undefined
860 USE quadric
861 USE quadric
862 USE sendrecv
863
864 IMPLICIT NONE
865
866 DOUBLE PRECISION x1,x2,x3
867 DOUBLE PRECISION f
868 INTEGER :: Q,Q_ID,BCID
869 LOGICAL :: CLIP_FLAG
870 CHARACTER (LEN = 7) :: METHOD
871 CHARACTER(LEN=9) :: GR
872
873 INTEGER :: GROUP,GS,P
874
875 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F_G
876
877 ALLOCATE(F_G(N_GROUP,0:DIM_QUADRIC))
878
879 SELECT CASE(METHOD)
880
881 CASE('QUADRIC')
882
883 F_G = - UNDEFINED
884
885 DO GROUP = 1, N_GROUP
886 GS = GROUP_SIZE(GROUP)
887 GR = TRIM(GROUP_RELATION(GROUP))
888
889 DO P = 1 , GS
890 Q_ID = GROUP_Q(GROUP,P)
891 CALL GET_F_QUADRIC(x1,x2,x3,Q_ID,F_G(GROUP,P),CLIP_FLAG)
892 ENDDO
893 IF(GR == 'AND') THEN
894 F_G(GROUP,0) = MAXVAL(F_G(GROUP,1:GS))
895 ELSEIF(GR == 'OR') THEN
896 F_G(GROUP,0) = MINVAL(F_G(GROUP,1:GS))
897 ELSEIF(GR == 'PIECEWISE') THEN
898 CALL REASSSIGN_QUADRIC(x1,x2,x3,GROUP,Q_ID)
899
900 CALL GET_F_QUADRIC(x1,x2,x3,Q_ID,F_G(GROUP,0),CLIP_FLAG)
901
902 ENDIF
903
904 ENDDO
905
906 f = F_G(1,0)
907
908 DO GROUP = 2, N_GROUP
909
910 GR = TRIM(RELATION_WITH_PREVIOUS(GROUP))
911
912 IF(GR =='AND') THEN
913 f = DMAX1(f,F_G(GROUP,0))
914 ELSEIF(GR =='OR') THEN
915 f = DMIN1(f,F_G(GROUP,0))
916 ENDIF
917
918 ENDDO
919
920
921 CASE('POLYGON')
922
923 CALL EVAL_POLY_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
924
925 CASE('USR_DEF')
926
927 CALL EVAL_USR_FCT(x1,x2,x3,Q,f,CLIP_FLAG)
928
929
930
931
932
933 CASE DEFAULT
934
935 WRITE(*,*)'ERROR IN SUBROUTINE EVAL_F.'
936 WRITE(*,*)'UNKNOWN METHOD:',METHOD
937 WRITE(*,*)'ACCEPTABLE METHODS:'
938 WRITE(*,*)'QUADRIC'
939 WRITE(*,*)'POLYGON'
940 WRITE(*,*)'USR_DEF'
941
942 CALL MFIX_EXIT(myPE)
943
944 END SELECT
945
946 DEALLOCATE(F_G)
947
948 RETURN
949 END SUBROUTINE EVAL_F
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965 SUBROUTINE INTERSECT_LINE(METHOD,xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
966
967 USE param
968 USE param1
969 USE parallel
970 USE constant
971 USE run
972 USE toleranc
973 USE geometry
974 USE indices
975 USE compar
976 USE sendrecv
977 USE quadric
978
979 IMPLICIT NONE
980 DOUBLE PRECISION:: x1,y1,z1,x2,y2,z2,x3,y3,z3
981 DOUBLE PRECISION:: xa,ya,za,xb,yb,zb,xc,yc,zc
982 INTEGER :: Q_ID,niter
983 DOUBLE PRECISION :: x_intersection
984 DOUBLE PRECISION :: f1,f2,f3,fa,fb
985 DOUBLE PRECISION :: t1,t2,t3
986 LOGICAL :: CLIP_FLAG,CLIP_FLAG1,CLIP_FLAG2,CLIP_FLAG3,INTERSECT_FLAG
987 CHARACTER (LEN=7) ::METHOD
988
989 x1 = xa
990 = ya
991 z1 = za
992 t1 = ZERO
993
994 x2 = xb
995 y2 = yb
996 z2 = zb
997 t2 = ONE
998
999 CALL EVAL_F(METHOD,x1,y1,z1,Q_ID,f1,CLIP_FLAG1)
1000 CALL EVAL_F(METHOD,x2,y2,z2,Q_ID,f2,CLIP_FLAG2)
1001
1002
1003
1004
1005
1006
1007 = 1
1008
1009 CLIP_FLAG = (CLIP_FLAG1).AND.(CLIP_FLAG2)
1010
1011 if(DABS(f1)<TOL_F) then
1012 = UNDEFINED
1013 yc = UNDEFINED
1014 zc = UNDEFINED
1015 INTERSECT_FLAG = .FALSE.
1016 elseif(DABS(f2)<TOL_F) then
1017 = UNDEFINED
1018 yc = UNDEFINED
1019 zc = UNDEFINED
1020 INTERSECT_FLAG = .FALSE.
1021 elseif(f1*f2 < ZERO) then
1022 niter = 0
1023 f3 = 2.0d0*TOL_F
1024 do while ( (abs(f3) > TOL_F) .AND. (niter<ITERMAX_INT) )
1025
1026 t3 = t1 - f1*(t2-t1)/(f2-f1)
1027
1028 x3 = x1 + t3 * (x2 - x1)
1029 y3 = y1 + t3 * (y2 - y1)
1030 z3 = z1 + t3 * (z2 - z1)
1031
1032 CALL EVAL_F(METHOD,x3,y3,z3,Q_ID,f3,CLIP_FLAG3)
1033 if(f1*f3<0) then
1034 t2 = t3
1035 f2 = f3
1036 else
1037 t1 = t3
1038 f1 = f3
1039 endif
1040 niter = niter + 1
1041
1042 end do
1043 if (niter < ITERMAX_INT) then
1044 xc = x3
1045 yc = y3
1046 zc = z3
1047 INTERSECT_FLAG = .TRUE.
1048 else
1049 WRITE(*,*)' Subroutine intersect_line:'
1050 WRITE(*,*) 'Unable to find the intersection of quadric:',Q_ID
1051 WRITE(*,1000)'between (x1,y1,z1)= ', xa,ya,za
1052 WRITE(*,1000)' and (x2,y2,z2)= ', xb,yb,zb
1053 CALL EVAL_F(METHOD,xa,ya,za,Q_ID,fa,CLIP_FLAG1)
1054 CALL EVAL_F(METHOD,xb,yb,zb,Q_ID,fb,CLIP_FLAG1)
1055 WRITE(*,1000)'f(x1,y1,z1) = ', fa
1056 WRITE(*,1000)'f(x2,y2,z2) = ', fb
1057 WRITE(*,1000)'Current Location (x3,y3,z3)= ', x3,y3,z3
1058 WRITE(*,1000)'Current value of abs(f) = ', DABS(f3)
1059 WRITE(*,1000)'Tolerance = ', TOL_F
1060 WRITE(*,*)'Maximum number of iterations = ', ITERMAX_INT
1061 WRITE(*,*) 'Please increase the intersection tolerance, '
1062 WRITE(*,*) 'or the maximum number of iterations, and try again.'
1063 WRITE(*,*) 'MFiX will exit now.'
1064 CALL MFIX_EXIT(myPE)
1065 x_intersection = UNDEFINED
1066 INTERSECT_FLAG = .FALSE.
1067
1068 endif
1069 else
1070 xc = UNDEFINED
1071 yc = UNDEFINED
1072 zc = UNDEFINED
1073 INTERSECT_FLAG = .FALSE.
1074 endif
1075
1076 1000 FORMAT(A,3(2X,G12.5))
1077
1078
1079 RETURN
1080
1081 END SUBROUTINE INTERSECT_LINE
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096 SUBROUTINE INTERSECT(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
1097
1098 USE param
1099 USE param1
1100 USE parallel
1101 USE constant
1102 USE run
1103 USE toleranc
1104 USE geometry
1105 USE indices
1106 USE compar
1107 USE sendrecv
1108 USE quadric
1109 USE cutcell
1110 USE polygon
1111 USE functions
1112
1113 IMPLICIT NONE
1114 CHARACTER (LEN=*) :: TYPE_OF_CELL
1115 INTEGER :: IJK,I,J,K,Q_ID,N_int_x,N_int_y,N_int_z,N_USR
1116 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb,xc,yc,zc
1117 DOUBLE PRECISION :: Xi,Yi,Zi,Xc_backup,Yc_backup,Zc_backup
1118 LOGICAL :: INTERSECT_FLAG
1119
1120 Xi = UNDEFINED
1121 Yi = UNDEFINED
1122 Zi = UNDEFINED
1123
1124
1125
1126
1127
1128 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1129
1130
1131
1132
1133 = X_NODE(7)
1134 ya = Y_NODE(7)
1135 za = Z_NODE(7)
1136
1137 xb = X_NODE(8)
1138 yb = Y_NODE(8)
1139 zb = Z_NODE(8)
1140
1141
1142 N_int_x = 0
1143 INTERSECT_X(IJK) = .FALSE.
1144 Q_ID = 1
1145 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
1146 IF ((INTERSECT_FLAG).AND.(xc/=Xi)) THEN
1147 N_int_x = N_int_x + 1
1148 INTERSECT_X(IJK) = .TRUE.
1149 xc_backup = Xi
1150 Xi = xc
1151 ENDIF
1152
1153 IF(N_int_x /= 1) THEN
1154 Xi = UNDEFINED
1155 INTERSECT_X(IJK) = .FALSE.
1156 ENDIF
1157
1158
1159 DO Q_ID = 1, N_POLYGON
1160 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_X(IJK),xc,yc,zc)
1161 IF(INTERSECT_X(IJK)) Xi = xc
1162 ENDDO
1163
1164 DO N_USR= 1, N_USR_DEF
1165 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_X(IJK),xc,yc,zc)
1166 IF(INTERSECT_X(IJK)) Xi = xc
1167 ENDDO
1168
1169
1170
1171
1172
1173
1174 IF(TYPE_OF_CELL=='U_MOMENTUM') THEN
1175 IF(SNAP(IJK)) THEN
1176 INTERSECT_X(IJK) = .TRUE.
1177 I = I_OF(IJK)
1178 Xi = XG_E(I)
1179 ENDIF
1180 ENDIF
1181
1182
1183
1184
1185
1186
1187 = X_NODE(6)
1188 ya = Y_NODE(6)
1189 za = Z_NODE(6)
1190
1191 N_int_y = 0
1192 INTERSECT_Y(IJK) = .FALSE.
1193 Q_ID = 1
1194 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
1195 IF ((INTERSECT_FLAG).AND.(yc/=Yi)) THEN
1196 N_int_y = N_int_y + 1
1197 INTERSECT_Y(IJK) = .TRUE.
1198 yc_backup = Yi
1199 Yi = yc
1200 ENDIF
1201
1202 IF(N_int_y /= 1) THEN
1203 Yi = UNDEFINED
1204 INTERSECT_Y(IJK) = .FALSE.
1205 ENDIF
1206
1207 DO Q_ID = 1, N_POLYGON
1208 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Y(IJK),xc,yc,zc)
1209 IF(INTERSECT_Y(IJK)) Yi = yc
1210 ENDDO
1211
1212 DO N_USR= 1, N_USR_DEF
1213 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Y(IJK),xc,yc,zc)
1214 IF(INTERSECT_Y(IJK)) Yi = yc
1215 ENDDO
1216
1217
1218
1219
1220
1221
1222 IF(TYPE_OF_CELL=='V_MOMENTUM') THEN
1223 IF(SNAP(IJK)) THEN
1224 INTERSECT_Y(IJK) = .TRUE.
1225 J = J_OF(IJK)
1226 Yi = YG_N(J)
1227 ENDIF
1228 ENDIF
1229
1230 IF(DO_K) THEN
1231
1232
1233
1234 = X_NODE(4)
1235 ya = Y_NODE(4)
1236 za = Z_NODE(4)
1237
1238 N_int_z = 0
1239 INTERSECT_Z(IJK) = .FALSE.
1240 Q_ID = 1
1241 CALL INTERSECT_LINE('QUADRIC',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_FLAG,xc,yc,zc)
1242 IF ((INTERSECT_FLAG).AND.(zc/=Zi)) THEN
1243 N_int_z = N_int_z + 1
1244 INTERSECT_Z(IJK) = .TRUE.
1245 zc_backup = Zi
1246 Zi = zc
1247 ENDIF
1248
1249 IF(N_int_z /= 1) THEN
1250 Zi = UNDEFINED
1251 INTERSECT_Z(IJK) = .FALSE.
1252 ENDIF
1253
1254 DO Q_ID = 1, N_POLYGON
1255 CALL INTERSECT_LINE('POLYGON',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Z(IJK),xc,yc,zc)
1256 IF(INTERSECT_Z(IJK)) Zi = zc
1257 ENDDO
1258
1259 DO N_USR= 1, N_USR_DEF
1260 CALL INTERSECT_LINE('USR_DEF',xa,ya,za,xb,yb,zb,Q_ID,INTERSECT_Z(IJK),xc,yc,zc)
1261 IF(INTERSECT_Z(IJK)) Zi = zc
1262 ENDDO
1263
1264
1265
1266
1267
1268
1269 IF(TYPE_OF_CELL=='W_MOMENTUM') THEN
1270 IF(SNAP(IJK)) THEN
1271 INTERSECT_Z(IJK) = .TRUE.
1272 K = K_OF(IJK)
1273 Zi = ZG_T(K)
1274 ENDIF
1275 ENDIF
1276
1277 ENDIF
1278
1279
1280 IF(INTERSECT_X(IJK)) THEN
1281 POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
1282 POTENTIAL_CUT_CELL_AT(EAST_OF(IJK)) = .TRUE.
1283 ENDIF
1284
1285 IF(INTERSECT_Y(IJK)) THEN
1286 POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
1287 POTENTIAL_CUT_CELL_AT(NORTH_OF(IJK)) = .TRUE.
1288 ENDIF
1289
1290
1291 IF(INTERSECT_Z(IJK)) THEN
1292 POTENTIAL_CUT_CELL_AT(IJK) = .TRUE.
1293 POTENTIAL_CUT_CELL_AT(TOP_OF(IJK)) = .TRUE.
1294 ENDIF
1295
1296
1297 RETURN
1298
1299 END SUBROUTINE INTERSECT
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316 SUBROUTINE CLEAN_INTERSECT(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
1317
1318 USE param
1319 USE param1
1320 USE parallel
1321 USE constant
1322 USE run
1323 USE toleranc
1324 USE geometry
1325 USE indices
1326 USE compar
1327 USE sendrecv
1328 USE quadric
1329 USE cutcell
1330 USE polygon
1331 USE STL
1332 USE functions
1333
1334 IMPLICIT NONE
1335 CHARACTER (LEN=*) :: TYPE_OF_CELL
1336 INTEGER :: IJK,I,J,K,IM,JM,KM,IP,JP,KP
1337 INTEGER :: BCID
1338 INTEGER :: IMJK,IPJK,IJMK,IJPK,IJKM,IJKP,IMJPK,IMJKP,IPJMK,IJMKP,IPJKM,IJPKM
1339 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb
1340 DOUBLE PRECISION :: Xi,Yi,Zi
1341 DOUBLE PRECISION :: DFC,DFC_MAX,F4,F6,F7,F8
1342 LOGICAL :: CLIP_FLAG,CAD,F_TEST
1343
1344
1345
1346
1347
1348
1349 = USE_MSH.OR.USE_STL
1350 F_TEST = .TRUE.
1351
1352
1353
1354
1355
1356 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1357
1358 I = I_OF(IJK)
1359 J = J_OF(IJK)
1360 K = K_OF(IJK)
1361
1362 IM = I - 1
1363 JM = J - 1
1364 KM = K - 1
1365
1366 IP = I + 1
1367 JP = J + 1
1368 KP = K + 1
1369
1370 IMJK = FUNIJK(IM,J,K)
1371 IPJK = FUNIJK(IP,J,K)
1372 IJMK = FUNIJK(I,JM,K)
1373 IJPK = FUNIJK(I,JP,K)
1374 IJKM = FUNIJK(I,J,KM)
1375 IJKP = FUNIJK(I,J,KP)
1376
1377 IMJPK = FUNIJK(IM,JP,K)
1378 IMJKP = FUNIJK(IM,J,KP)
1379
1380 IPJMK = FUNIJK(IP,JM,K)
1381 IJMKP = FUNIJK(I,JM,KP)
1382
1383 IPJKM = FUNIJK(IP,J,KM)
1384 IJPKM = FUNIJK(I,JP,KM)
1385
1386
1387 IF(IMJK<1.OR.IMJK>DIMENSION_3) IMJK = IJK
1388 IF(IPJK<1.OR.IPJK>DIMENSION_3) IPJK = IJK
1389 IF(IJMK<1.OR.IJMK>DIMENSION_3) IJMK = IJK
1390 IF(IJPK<1.OR.IJPK>DIMENSION_3) IJPK = IJK
1391 IF(IJKM<1.OR.IJKM>DIMENSION_3) IJKM = IJK
1392 IF(IJKP<1.OR.IJKP>DIMENSION_3) IJKP = IJK
1393
1394 IF(IMJPK<1.OR.IMJPK>DIMENSION_3) IMJPK = IJK
1395 IF(IMJKP<1.OR.IMJKP>DIMENSION_3) IMJKP = IJK
1396
1397 IF(IPJMK<1.OR.IPJMK>DIMENSION_3) IPJMK = IJK
1398 IF(IJMKP<1.OR.IJMKP>DIMENSION_3) IJMKP = IJK
1399
1400 IF(IPJKM<1.OR.IPJKM>DIMENSION_3) IPJKM = IJK
1401 IF(IJPKM<1.OR.IJPKM>DIMENSION_3) IJPKM = IJK
1402
1403
1404
1405
1406
1407 = X_NODE(7)
1408 ya = Y_NODE(7)
1409 za = Z_NODE(7)
1410
1411 xb = X_NODE(8)
1412 yb = Y_NODE(8)
1413 zb = Z_NODE(8)
1414
1415 DFC_MAX = TOL_SNAP(1) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1416
1417 IF(INTERSECT_X(IJK)) THEN
1418
1419 DFC = DABS(Xi-xa)
1420
1421 IF(CAD) F_TEST = (F_AT(IMJK)/=ZERO)
1422 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1423 IF(PRINT_WARNINGS) THEN
1424 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 7'
1425 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1426 ENDIF
1427
1428 INTERSECT_X(IJK) = .FALSE.
1429 IF(I>=IMIN1) THEN
1430 INTERSECT_X(IMJK) = .FALSE.
1431 INTERSECT_Y(IMJK) = .FALSE.
1432 INTERSECT_Y(IMJPK) = .FALSE.
1433 IF(DO_K) INTERSECT_Z(IMJK) = .FALSE.
1434 IF(DO_K) INTERSECT_Z(IMJKP) = .FALSE.
1435
1436 SNAP(IMJK) = .TRUE.
1437 ENDIF
1438
1439 ENDIF
1440
1441
1442 DFC = DABS(Xi-xb)
1443
1444 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1445 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1446 IF(PRINT_WARNINGS) THEN
1447 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 8'
1448 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1449 ENDIF
1450
1451
1452 INTERSECT_X(IJK) = .FALSE.
1453 INTERSECT_Y(IJK) = .FALSE.
1454 IF(DO_K) INTERSECT_Z(IJK) = .FALSE.
1455 SNAP(IJK) = .TRUE.
1456
1457 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
1458 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
1459 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
1460
1461
1462 ENDIF
1463
1464 IF(USE_STL.OR.USE_MSH) THEN
1465 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,7,F7,CLIP_FLAG,BCID)
1466 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1467 IF(F7*F8>TOL_STL**2) INTERSECT_X(IJK) = .FALSE.
1468 ENDIF
1469
1470 ENDIF
1471
1472
1473
1474
1475
1476
1477
1478
1479 = X_NODE(6)
1480 ya = Y_NODE(6)
1481 za = Z_NODE(6)
1482
1483 DFC_MAX = TOL_SNAP(2) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1484
1485
1486 IF(INTERSECT_Y(IJK)) THEN
1487
1488 DFC = DABS(Yi-ya)
1489
1490 IF(CAD) F_TEST = (F_AT(IJMK)/=ZERO)
1491 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1492 IF(PRINT_WARNINGS) THEN
1493 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 6'
1494 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1495 ENDIF
1496
1497 INTERSECT_Y(IJK) = .FALSE.
1498 IF(J>=JMIN1) THEN
1499 INTERSECT_X(IJMK) = .FALSE.
1500 INTERSECT_X(IPJMK) = .FALSE.
1501 INTERSECT_Y(IJMK) = .FALSE.
1502 IF(DO_K) INTERSECT_Z(IJMK) = .FALSE.
1503 IF(DO_K) INTERSECT_Z(IJMKP) = .FALSE.
1504
1505 SNAP(IJMK) = .TRUE.
1506 ENDIF
1507
1508 ENDIF
1509
1510
1511 DFC = DABS(Yi-yb)
1512
1513 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1514 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1515 IF(PRINT_WARNINGS) THEN
1516 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 8'
1517 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1518 ENDIF
1519
1520 INTERSECT_X(IJK) = .FALSE.
1521 INTERSECT_Y(IJK) = .FALSE.
1522 IF(DO_K) INTERSECT_Z(IJK) = .FALSE.
1523 SNAP(IJK) = .TRUE.
1524
1525 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
1526 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
1527 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
1528
1529
1530 ENDIF
1531
1532
1533 IF(USE_STL.OR.USE_MSH) THEN
1534 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,6,F6,CLIP_FLAG,BCID)
1535 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1536 IF(F6*F8>TOL_STL**2) INTERSECT_Y(IJK) = .FALSE.
1537 ENDIF
1538
1539 ENDIF
1540
1541
1542 IF(DO_K) THEN
1543
1544
1545
1546
1547 = X_NODE(4)
1548 ya = Y_NODE(4)
1549 za = Z_NODE(4)
1550
1551 DFC_MAX = TOL_SNAP(3) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1552
1553 IF(INTERSECT_Z(IJK)) THEN
1554
1555 DFC = DABS(Zi-Za)
1556
1557 IF(CAD) F_TEST = (F_AT(IJKM)/=ZERO)
1558 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1559 IF(PRINT_WARNINGS) THEN
1560 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 4'
1561 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1562 ENDIF
1563
1564 INTERSECT_Z(IJK) = .FALSE.
1565
1566 IF(K>=KMIN1) THEN
1567 INTERSECT_X(IJKM) = .FALSE.
1568 INTERSECT_X(IPJKM) = .FALSE.
1569 INTERSECT_Y(IJKM) = .FALSE.
1570 INTERSECT_Y(IJPKM) = .FALSE.
1571 INTERSECT_Z(IJKM) = .FALSE.
1572
1573 SNAP(IJKM) = .TRUE.
1574 ENDIF
1575
1576 ENDIF
1577
1578
1579 DFC = DABS(Zi-Zb)
1580
1581 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1582 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1583 IF(PRINT_WARNINGS) THEN
1584 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 8'
1585 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1586 ENDIF
1587
1588 INTERSECT_X(IJK) = .FALSE.
1589 INTERSECT_Y(IJK) = .FALSE.
1590 INTERSECT_Z(IJK) = .FALSE.
1591 SNAP(IJK) = .TRUE.
1592
1593
1594
1595
1596 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
1597 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
1598 IF(K<=KMAX1) INTERSECT_Z(IJKP) = .FALSE.
1599
1600
1601 ENDIF
1602
1603 IF(USE_STL.OR.USE_MSH) THEN
1604 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,4,F4,CLIP_FLAG,BCID)
1605 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1606 IF(F4*F8>TOL_STL**2) INTERSECT_Z(IJK) = .FALSE.
1607 ENDIF
1608
1609 ENDIF
1610
1611 ENDIF
1612
1613
1614 RETURN
1615
1616 END SUBROUTINE CLEAN_INTERSECT
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633 SUBROUTINE CLEAN_INTERSECT_SCALAR
1634
1635 USE param
1636 USE param1
1637 USE parallel
1638 USE constant
1639 USE run
1640 USE toleranc
1641 USE geometry
1642 USE indices
1643 USE compar
1644 USE sendrecv
1645 USE quadric
1646 USE cutcell
1647 USE polygon
1648 USE STL
1649 USE functions
1650
1651 IMPLICIT NONE
1652 INTEGER :: IJK
1653
1654
1655
1656 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1657 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1658 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1659 call send_recv(Xn_int,2)
1660 call send_recv(Ye_int,2)
1661 call send_recv(Zt_int,2)
1662
1663 DO IJK = IJKSTART3, IJKEND3
1664 CALL SET_SNAP_FLAG(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
1665 END DO
1666
1667 call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1668 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1669 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1670 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1671 call send_recv(Xn_int,2)
1672 call send_recv(Ye_int,2)
1673 call send_recv(Zt_int,2)
1674
1675 DO IJK = IJKSTART3, IJKEND3
1676 CALL REMOVE_INTERSECT_FLAG(IJK)
1677 END DO
1678
1679 call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
1680 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
1681 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
1682 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
1683 call send_recv(Xn_int,2)
1684 call send_recv(Ye_int,2)
1685 call send_recv(Zt_int,2)
1686
1687 RETURN
1688
1689 END SUBROUTINE CLEAN_INTERSECT_SCALAR
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705 SUBROUTINE SET_SNAP_FLAG(IJK,TYPE_OF_CELL,Xi,Yi,Zi)
1706
1707 USE param
1708 USE param1
1709 USE parallel
1710 USE constant
1711 USE run
1712 USE toleranc
1713 USE geometry
1714 USE indices
1715 USE compar
1716 USE sendrecv
1717 USE quadric
1718 USE cutcell
1719 USE polygon
1720 USE STL
1721 USE functions
1722
1723 IMPLICIT NONE
1724 CHARACTER (LEN=*) :: TYPE_OF_CELL
1725 INTEGER :: IJK,I,J,K,IM,JM,KM
1726 INTEGER :: BCID
1727 INTEGER :: IMJK,IJMK,IJKM
1728 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb
1729 DOUBLE PRECISION :: Xi,Yi,Zi
1730 DOUBLE PRECISION :: DFC,DFC_MAX,F4,F6,F7,F8
1731 LOGICAL :: CLIP_FLAG,CAD,F_TEST
1732
1733
1734
1735
1736
1737
1738 = USE_MSH.OR.USE_STL
1739 F_TEST = .TRUE.
1740
1741
1742
1743
1744
1745 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
1746
1747 I = I_OF(IJK)
1748 J = J_OF(IJK)
1749 K = K_OF(IJK)
1750
1751 IM = I - 1
1752 JM = J - 1
1753 KM = K - 1
1754
1755 IMJK = FUNIJK(IM,J,K)
1756 IJMK = FUNIJK(I,JM,K)
1757 IJKM = FUNIJK(I,J,KM)
1758
1759 IF(IMJK<1.OR.IMJK>DIMENSION_3) IMJK = IJK
1760 IF(IJMK<1.OR.IJMK>DIMENSION_3) IJMK = IJK
1761 IF(IJKM<1.OR.IJKM>DIMENSION_3) IJKM = IJK
1762
1763
1764
1765
1766
1767 = X_NODE(7)
1768 ya = Y_NODE(7)
1769 za = Z_NODE(7)
1770
1771 xb = X_NODE(8)
1772 yb = Y_NODE(8)
1773 zb = Z_NODE(8)
1774
1775
1776 DFC_MAX = TOL_SNAP(1) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1777
1778 IF(INTERSECT_X(IJK)) THEN
1779
1780 DFC = DABS(Xi-xa)
1781
1782 IF(CAD) F_TEST = (F_AT(IMJK)/=ZERO)
1783 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1784 IF(PRINT_WARNINGS) THEN
1785 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 7'
1786 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1787 ENDIF
1788
1789 IF(I>=IMIN1) SNAP(IMJK) = .TRUE.
1790
1791 ENDIF
1792
1793
1794 DFC = DABS(Xi-xb)
1795
1796 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1797 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1798 IF(PRINT_WARNINGS) THEN
1799 WRITE(*,*)'MERGING X-INTERSECTION ALONG EDGE 7 ONTO NODE 8'
1800 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1801 ENDIF
1802
1803 SNAP(IJK) = .TRUE.
1804
1805 ENDIF
1806
1807 IF(USE_STL.OR.USE_MSH) THEN
1808 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,7,F7,CLIP_FLAG,BCID)
1809 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1810 IF(F7*F8>TOL_STL**2) INTERSECT_X(IJK) = .FALSE.
1811 ENDIF
1812
1813 ENDIF
1814
1815
1816
1817
1818
1819
1820 = X_NODE(6)
1821 ya = Y_NODE(6)
1822 za = Z_NODE(6)
1823
1824 DFC_MAX = TOL_SNAP(2) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1825
1826
1827 IF(INTERSECT_Y(IJK)) THEN
1828
1829 DFC = DABS(Yi-ya)
1830
1831 IF(CAD) F_TEST = (F_AT(IJMK)/=ZERO)
1832 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1833 IF(PRINT_WARNINGS) THEN
1834 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 6'
1835 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1836 ENDIF
1837
1838 IF(J>=JMIN1) SNAP(IJMK) = .TRUE.
1839
1840 ENDIF
1841
1842
1843 DFC = DABS(Yi-yb)
1844
1845 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1846 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1847 IF(PRINT_WARNINGS) THEN
1848 WRITE(*,*)'MERGING Y-INTERSECTION ALONG EDGE 6 ONTO NODE 8'
1849 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1850 ENDIF
1851
1852 SNAP(IJK) = .TRUE.
1853
1854
1855 ENDIF
1856
1857
1858 IF(USE_STL.OR.USE_MSH) THEN
1859 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,6,F6,CLIP_FLAG,BCID)
1860 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1861 IF(F6*F8>TOL_STL**2) INTERSECT_Y(IJK) = .FALSE.
1862 ENDIF
1863
1864 ENDIF
1865
1866
1867 IF(DO_K) THEN
1868
1869
1870
1871
1872 = X_NODE(4)
1873 ya = Y_NODE(4)
1874 za = Z_NODE(4)
1875
1876 DFC_MAX = TOL_SNAP(3) * DSQRT((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
1877
1878 IF(INTERSECT_Z(IJK)) THEN
1879
1880 DFC = DABS(Zi-Za)
1881
1882 IF(CAD) F_TEST = (F_AT(IJKM)/=ZERO)
1883 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1884 IF(PRINT_WARNINGS) THEN
1885 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 4'
1886 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1887 ENDIF
1888
1889 IF(K>=KMIN1) SNAP(IJKM) = .TRUE.
1890
1891 ENDIF
1892
1893 DFC = DABS(Zi-Zb)
1894
1895 IF(CAD) F_TEST = (F_AT(IJK)/=ZERO)
1896 IF(DFC < DFC_MAX.AND.F_TEST) THEN
1897 IF(PRINT_WARNINGS) THEN
1898 WRITE(*,*)'MERGING Z-INTERSECTION ALONG EDGE 11 ONTO NODE 8'
1899 WRITE(*,*)'AT IJK,I,J,K=',IJK,I,J,K
1900 ENDIF
1901
1902 SNAP(IJK) = .TRUE.
1903
1904 ENDIF
1905
1906
1907 IF(USE_STL.OR.USE_MSH) THEN
1908 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,4,F4,CLIP_FLAG,BCID)
1909 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,8,F8,CLIP_FLAG,BCID)
1910 IF(F4*F8>TOL_STL**2) INTERSECT_Z(IJK) = .FALSE.
1911 ENDIF
1912
1913 ENDIF
1914
1915 ENDIF
1916
1917
1918 RETURN
1919
1920 END SUBROUTINE SET_SNAP_FLAG
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936 SUBROUTINE REMOVE_INTERSECT_FLAG(IJK)
1937
1938 USE param
1939 USE param1
1940 USE parallel
1941 USE constant
1942 USE run
1943 USE toleranc
1944 USE geometry
1945 USE indices
1946 USE compar
1947 USE sendrecv
1948 USE quadric
1949 USE cutcell
1950 USE polygon
1951 USE STL
1952 USE functions
1953
1954 IMPLICIT NONE
1955 INTEGER :: IJK,I,J,K,IP,JP,KP
1956 INTEGER :: IPJK,IJPK,IJKP
1957
1958
1959 I = I_OF(IJK)
1960 J = J_OF(IJK)
1961 K = K_OF(IJK)
1962
1963 IP = I + 1
1964 JP = J + 1
1965 KP = K + 1
1966
1967 IPJK = FUNIJK(IP,J,K)
1968 IJPK = FUNIJK(I,JP,K)
1969 IJKP = FUNIJK(I,J,KP)
1970
1971 IF(IPJK<1.OR.IPJK>DIMENSION_3) IPJK = IJK
1972 IF(IJPK<1.OR.IJPK>DIMENSION_3) IJPK = IJK
1973 IF(IJKP<1.OR.IJKP>DIMENSION_3) IJKP = IJK
1974
1975 IF(SNAP(IJK)) THEN
1976 INTERSECT_X(IJK) = .FALSE.
1977 INTERSECT_Y(IJK) = .FALSE.
1978 INTERSECT_Z(IJK) = .FALSE.
1979 IF(I<=IMAX1) INTERSECT_X(IPJK) = .FALSE.
1980 IF(J<=JMAX1) INTERSECT_Y(IJPK) = .FALSE.
1981 IF(DO_K.AND.(K<=KMAX1)) INTERSECT_Z(IJKP) = .FALSE.
1982 ENDIF
1983
1984
1985 RETURN
1986
1987 END SUBROUTINE REMOVE_INTERSECT_FLAG
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003 SUBROUTINE OPEN_CUT_CELL_FILES
2004
2005 USE cutcell
2006 USE compar
2007
2008 IMPLICIT NONE
2009
2010 IF(MyPE == PE_IO) THEN
2011 OPEN(CONVERT='BIG_ENDIAN',UNIT = UNIT_CUT_CELL_LOG, FILE = 'CUT_CELL.LOG')
2012 ENDIF
2013
2014 RETURN
2015
2016
2017 END SUBROUTINE OPEN_CUT_CELL_FILES
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032 SUBROUTINE CLOSE_CUT_CELL_FILES
2033
2034 USE compar
2035 USE cutcell
2036
2037 IMPLICIT NONE
2038
2039 IF(MyPE == PE_IO) THEN
2040 CLOSE(UNIT_CUT_CELL_LOG)
2041 ENDIF
2042
2043 RETURN
2044
2045
2046 END SUBROUTINE CLOSE_CUT_CELL_FILES
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063 SUBROUTINE CAD_INTERSECT(TYPE_OF_CELL,Xint,Yint,Zint)
2064
2065 USE param
2066 USE param1
2067 USE parallel
2068 USE constant
2069 USE run
2070 USE toleranc
2071 USE geometry
2072 USE indices
2073 USE compar
2074 USE sendrecv
2075 USE quadric
2076 USE cutcell
2077 USE polygon
2078 USE stl
2079
2080 USE mpi_utility
2081 USE functions
2082
2083 IMPLICIT NONE
2084 CHARACTER (LEN=*) :: TYPE_OF_CELL
2085 INTEGER :: IJK,I,J,K
2086 INTEGER :: IM,IP,JM,JP,KM,KP,IMJK,IPJK,IJMK,IJPK,IJKM,IJKP
2087 INTEGER :: IJPKP,IPJKP,IPJPK
2088 DOUBLE PRECISION :: xa,ya,za,xb,yb,zb,xc,yc,zc
2089 LOGICAL :: INTERSECT_FLAG,INSIDE_FACET_a,INSIDE_FACET_b
2090
2091 DOUBLE PRECISION :: X1,X2,Y1,Y2,Z1,Z2
2092
2093 DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: Xint,Yint,Zint
2094
2095 INTEGER :: NN,I1,I2,J1,J2,K1,K2
2096
2097 DOUBLE PRECISION :: X_OFFSET, Y_OFFSET, Z_OFFSET
2098
2099 DOUBLE PRECISION, DIMENSION(3) :: N4,N6,N7,N8
2100 DOUBLE PRECISION :: CURRENT_F
2101
2102 INTEGER :: N_UNDEFINED, NTOTAL_UNDEFINED,N_PROP
2103 INTEGER, PARAMETER :: N_PROPMAX=1000
2104 LOGICAL:: F_FOUND
2105
2106
2107
2108 = .FALSE.
2109 INTERSECT_Y = .FALSE.
2110 INTERSECT_Z = .FALSE.
2111
2112 Xint = UNDEFINED
2113 Yint = UNDEFINED
2114 Zint = UNDEFINED
2115
2116 F_AT = UNDEFINED
2117
2118 SELECT CASE (TYPE_OF_CELL)
2119 CASE('SCALAR')
2120
2121
2122 X_OFFSET = ZERO
2123 Y_OFFSET = ZERO
2124 Z_OFFSET = ZERO
2125
2126 CASE('U_MOMENTUM')
2127
2128 X_OFFSET = HALF
2129 Y_OFFSET = ZERO
2130 Z_OFFSET = ZERO
2131
2132 CASE('V_MOMENTUM')
2133
2134 X_OFFSET = ZERO
2135 Y_OFFSET = HALF
2136 Z_OFFSET = ZERO
2137
2138 CASE('W_MOMENTUM')
2139
2140 X_OFFSET = ZERO
2141 Y_OFFSET = ZERO
2142 Z_OFFSET = HALF
2143
2144
2145 CASE DEFAULT
2146 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
2147 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
2148 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
2149 WRITE(*,*)'SCALAR'
2150 WRITE(*,*)'U_MOMENTUM'
2151 WRITE(*,*)'V_MOMENTUM'
2152 WRITE(*,*)'W_MOMENTUM'
2153 CALL MFIX_EXIT(myPE)
2154 END SELECT
2155
2156
2157 DO NN = 1,N_FACETS
2158
2159
2160 X1 = MINVAL(VERTEX(1:3,1,NN))
2161 X2 = MAXVAL(VERTEX(1:3,1,NN))
2162 Y1 = MINVAL(VERTEX(1:3,2,NN))
2163 Y2 = MAXVAL(VERTEX(1:3,2,NN))
2164 Z1 = MINVAL(VERTEX(1:3,3,NN))
2165 Z2 = MAXVAL(VERTEX(1:3,3,NN))
2166
2167
2168 I1 = IEND3
2169 I2 = ISTART3
2170
2171 IF(X2>=ZERO-DX(ISTART3).AND.X1<=XLENGTH+DX(IEND3)) THEN
2172 DO I = ISTART3, IEND3
2173 IP = I+1
2174 IF(XG_E(I)+X_OFFSET*DX(IP)>=X1-TOL_STL) THEN
2175 I1=I
2176 EXIT
2177 ENDIF
2178 ENDDO
2179
2180 DO I = IEND3, ISTART3,-1
2181 IP = I+1
2182 IF(XG_E(I)-DX(I)+X_OFFSET*DX(IP)<=X2+TOL_STL) THEN
2183 I2=I
2184 EXIT
2185 ENDIF
2186 ENDDO
2187 ENDIF
2188
2189
2190 J1 = JEND3
2191 J2 = JSTART3
2192
2193 IF(Y2>=ZERO-DY(JSTART3).AND.Y1<=YLENGTH+DY(JEND3)) THEN
2194 DO J = JSTART3, JEND3
2195 JP = J+1
2196 IF(YG_N(J)+Y_OFFSET*DY(JP)>=Y1-TOL_STL) THEN
2197 J1=J
2198 EXIT
2199 ENDIF
2200 ENDDO
2201
2202 DO J = JEND3, JSTART3,-1
2203 JP=J+1
2204 IF(YG_N(J)-DY(J)+Y_OFFSET*DY(JP)<=Y2+TOL_STL) THEN
2205 J2=J
2206 EXIT
2207 ENDIF
2208 ENDDO
2209 ENDIF
2210
2211 K1 = KEND3
2212 K2 = KSTART3
2213
2214 IF(Z2>=ZERO-DZ(KSTART3).AND.Z1<=ZLENGTH+DZ(KEND3)) THEN
2215 DO K = KSTART3, KEND3
2216 KP=K+1
2217
2218 IF(ZG_T(K)+Z_OFFSET*DZ(KP)>=Z1-TOL_STL) THEN
2219 K1=K
2220 EXIT
2221 ENDIF
2222 ENDDO
2223
2224 DO K = KEND3, KSTART3,-1
2225 KP = K+1
2226 IF(ZG_T(K)-DZ(K)+Z_OFFSET*DZ(KP)<=Z2+TOL_STL) THEN
2227 K2=K
2228 EXIT
2229 ENDIF
2230 ENDDO
2231 ENDIF
2232
2233
2234
2235
2236 DO K=K1,K2
2237 DO J=J1,J2
2238 DO I=I1,I2
2239
2240 IJK = FUNIJK(I,J,K)
2241
2242
2243 IM = MAX0(I - 1 ,ISTART3)
2244 JM = MAX0(J - 1 ,JSTART3)
2245 KM = MAX0(K - 1 ,KSTART3)
2246
2247 IP = MIN0(I + 1 ,IEND3)
2248 JP = MIN0(J + 1 ,JEND3)
2249 KP = MIN0(K + 1 ,KEND3)
2250
2251
2252 IMJK = FUNIJK(IM,J,K)
2253 IPJK = FUNIJK(IP,J,K)
2254 IJMK = FUNIJK(I,JM,K)
2255 IJPK = FUNIJK(I,JP,K)
2256 IJKM = FUNIJK(I,J,KM)
2257 IJKP = FUNIJK(I,J,KP)
2258
2259 IJPKP = FUNIJK(I,JP,KP)
2260 IPJKP = FUNIJK(IP,J,KP)
2261 IPJPK = FUNIJK(IP,JP,K)
2262
2263
2264
2265
2266
2267
2268 CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
2269
2270
2271
2272
2273 = X_NODE(7)
2274 ya = Y_NODE(7)
2275 za = Z_NODE(7)
2276
2277 xb = X_NODE(8)
2278 yb = Y_NODE(8)
2279 zb = Z_NODE(8)
2280
2281
2282
2283 CALL IS_POINT_INSIDE_FACET(xa,ya,za,NN,INSIDE_FACET_a)
2284
2285 IF(INSIDE_FACET_a) THEN
2286
2287 (IMJK) = ZERO
2288
2289 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2290
2291 ENDIF
2292
2293
2294 CALL IS_POINT_INSIDE_FACET(xb,yb,zb,NN,INSIDE_FACET_b)
2295
2296 IF(INSIDE_FACET_b) THEN
2297
2298 (IJK) = ZERO
2299
2300 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2301
2302 ENDIF
2303
2304
2305
2306
2307
2308 = .FALSE.
2309
2310 IF(.NOT.(INTERSECT_X(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
2311 CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,NN,INTERSECT_FLAG,xc,yc,zc)
2312 ENDIF
2313
2314 IF(INTERSECT_FLAG) THEN
2315 IF(INTERSECT_X(IJK)) THEN
2316 IF(DABS(Xint(IJK)-xc)>TOL_STL) THEN
2317
2318
2319 (IJK) = .FALSE.
2320 ENDIF
2321 ELSE
2322 INTERSECT_X(IJK) = .TRUE.
2323 Xint(IJK) = xc
2324
2325
2326
2327 (1) = xa-xc
2328 N7(2) = ya-yc
2329 N7(3) = za-zc
2330
2331 IF(DABS(F_AT(IMJK))>TOL_F) F_AT(IMJK) = -DOT_PRODUCT(N7,NORM_FACE(:,NN))
2332
2333 N8(1) = xb-xc
2334 N8(2) = yb-yc
2335 N8(3) = zb-zc
2336
2337 IF(DABS(F_AT(IJK))>TOL_F) F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,NN))
2338
2339
2340 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2341 IF(JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPK,NN)
2342 IF(KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJKP,NN)
2343 IF(JP<=J2.AND.KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPKP,NN)
2344 ENDIF
2345 ENDIF
2346
2347 IF(TYPE_OF_CELL=='U_MOMENTUM') THEN
2348 IF(SNAP(IJK)) THEN
2349 INTERSECT_X(IJK) = .TRUE.
2350 Xn_int(IJK) = XG_E(I)
2351 ENDIF
2352 ENDIF
2353
2354
2355
2356
2357 = X_NODE(6)
2358 ya = Y_NODE(6)
2359 za = Z_NODE(6)
2360
2361
2362
2363 CALL IS_POINT_INSIDE_FACET(xa,ya,za,NN,INSIDE_FACET_a)
2364
2365 IF(INSIDE_FACET_a) THEN
2366
2367 (IJMK) = ZERO
2368
2369 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2370
2371 ENDIF
2372
2373
2374
2375
2376
2377
2378
2379 = .FALSE.
2380
2381 IF(.NOT.(INTERSECT_Y(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
2382 CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,NN,INTERSECT_FLAG,xc,yc,zc)
2383 ENDIF
2384
2385
2386 IF(INTERSECT_FLAG) THEN
2387
2388 IF(INTERSECT_Y(IJK)) THEN
2389
2390 IF(DABS(Yint(IJK)-yc)>TOL_STL) THEN
2391
2392 INTERSECT_Y(IJK) = .FALSE.
2393
2394 ENDIF
2395
2396 ELSE
2397
2398
2399 INTERSECT_Y(IJK) = .TRUE.
2400 Yint(IJK) = yc
2401
2402
2403
2404 (1) = xa-xc
2405 N6(2) = ya-yc
2406 N6(3) = za-zc
2407
2408 IF(DABS(F_AT(IJMK))>TOL_F) F_AT(IJMK) = -DOT_PRODUCT(N6,NORM_FACE(:,NN))
2409
2410 N8(1) = xb-xc
2411 N8(2) = yb-yc
2412 N8(3) = zb-zc
2413
2414 IF(DABS(F_AT(IJK))>TOL_F) F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,NN))
2415
2416
2417 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2418 IF(IP<=I2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJK,NN)
2419 IF(KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJKP,NN)
2420 IF(IP<=I2.AND.KP<=K2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJKP,NN)
2421
2422 ENDIF
2423
2424 ENDIF
2425
2426 IF(TYPE_OF_CELL=='V_MOMENTUM') THEN
2427 IF(SNAP(IJK)) THEN
2428 INTERSECT_Y(IJK) = .TRUE.
2429 Ye_int(IJK) = YG_N(J)
2430 ENDIF
2431 ENDIF
2432
2433 IF(DO_K) THEN
2434
2435
2436
2437 = X_NODE(4)
2438 ya = Y_NODE(4)
2439 za = Z_NODE(4)
2440
2441
2442
2443 CALL IS_POINT_INSIDE_FACET(xa,ya,za,NN,INSIDE_FACET_a)
2444
2445 IF(INSIDE_FACET_a) THEN
2446
2447 (IJKM) = ZERO
2448
2449 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2450
2451 ENDIF
2452
2453
2454
2455
2456
2457
2458
2459 = .FALSE.
2460
2461 IF(.NOT.(INTERSECT_Z(IJK).OR.INSIDE_FACET_a.OR.INSIDE_FACET_b)) THEN
2462 CALL INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,NN,INTERSECT_FLAG,xc,yc,zc)
2463 ENDIF
2464
2465 IF(INTERSECT_FLAG) THEN
2466
2467 IF(INTERSECT_Z(IJK)) THEN
2468
2469 IF(DABS(Zint(IJK)-zc)>TOL_STL) THEN
2470
2471 INTERSECT_Z(IJK) = .FALSE.
2472
2473 ENDIF
2474
2475 ELSE
2476
2477
2478 INTERSECT_Z(IJK) = .TRUE.
2479 Zint(IJK) = zc
2480
2481
2482
2483
2484 (1) = xa-xc
2485 N4(2) = ya-yc
2486 N4(3) = za-zc
2487
2488 IF(DABS(F_AT(IJKM))>TOL_F) F_AT(IJKM) = -DOT_PRODUCT(N4,NORM_FACE(:,NN))
2489
2490 N8(1) = xb-xc
2491 N8(2) = yb-yc
2492 N8(3) = zb-zc
2493
2494 IF(DABS(F_AT(IJK))>TOL_F) F_AT(IJK) = -DOT_PRODUCT(N8,NORM_FACE(:,NN))
2495
2496
2497
2498 IF(TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJK,NN)
2499 IF(IP<=I2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJK,NN)
2500 IF(JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IJPK,NN)
2501 IF(IP<=I2.AND.JP<=J2.AND.TRIM(TYPE_OF_CELL).eq.'SCALAR') CALL ADD_FACET_AND_SET_BC_ID(IPJPK,NN)
2502
2503 ENDIF
2504
2505 ENDIF
2506
2507
2508 IF(TYPE_OF_CELL=='W_MOMENTUM') THEN
2509 IF(SNAP(IJK)) THEN
2510 INTERSECT_Z(IJK) = .TRUE.
2511 Zt_int(IJK) = ZG_T(K)
2512 ENDIF
2513 ENDIF
2514
2515 ENDIF
2516
2517
2518 ENDDO
2519 ENDDO
2520 ENDDO
2521
2522
2523 ENDDO
2524
2525 = UNDEFINED
2526
2527
2528
2529 DO IJK = IJKSTART3, IJKEND3
2530 IF(DABS(F_AT(IJK))<TOL_STL) THEN
2531 F_AT(IJK)=ZERO
2532 ENDIF
2533 END DO
2534
2535
2536
2537
2538
2539 call send_recv(F_AT,2)
2540 call send_recv(Xint,2)
2541 call send_recv(Yint,2)
2542 call send_recv(Zint,2)
2543 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
2544 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
2545 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
2546
2547
2548
2549
2550 DO IJK = IJKSTART3, IJKEND3
2551
2552
2553
2554
2555
2556 CALL CLEAN_INTERSECT(IJK,TYPE_OF_CELL,Xint(IJK),Yint(IJK),Zint(IJK))
2557
2558
2559
2560 END DO
2561
2562 call send_recv(F_AT,2)
2563 call SEND_RECEIVE_1D_LOGICAL(SNAP,2)
2564 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_X,2)
2565 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Y,2)
2566 call SEND_RECEIVE_1D_LOGICAL(INTERSECT_Z,2)
2567
2568 SELECT CASE (CAD_PROPAGATE_ORDER)
2569
2570 CASE (' ')
2571
2572
2573
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588 DO N_PROP=1,N_PROPMAX
2589
2590 DO IJK = IJKSTART3, IJKEND3
2591
2592
2593 = I_OF(IJK)
2594 J = J_OF(IJK)
2595 K = K_OF(IJK)
2596 IF(.NOT.IS_ON_myPE_plus1layer(I,J,K))cycle
2597
2598
2599 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2600
2601 IMJK = IM_OF(IJK)
2602 IF(F_AT(IMJK)==UNDEFINED.AND.F_AT(IMJK)/=ZERO) F_AT(IMJK)=F_AT(IJK)
2603
2604 IPJK = IP_OF(IJK)
2605 IF(F_AT(IPJK)==UNDEFINED.AND.F_AT(IPJK)/=ZERO) F_AT(IPJK)=F_AT(IJK)
2606
2607 IJMK = JM_OF(IJK)
2608 IF(F_AT(IJMK)==UNDEFINED.AND.F_AT(IJMK)/=ZERO) F_AT(IJMK)=F_AT(IJK)
2609
2610 IJPK = JP_OF(IJK)
2611 IF(F_AT(IJPK)==UNDEFINED.AND.F_AT(IJPK)/=ZERO) F_AT(IJPK)=F_AT(IJK)
2612
2613 IJKM = KM_OF(IJK)
2614 IF(F_AT(IJKM)==UNDEFINED.AND.F_AT(IJKM)/=ZERO) F_AT(IJKM)=F_AT(IJK)
2615
2616 IJKP = KP_OF(IJK)
2617 IF(F_AT(IJKP)==UNDEFINED.AND.F_AT(IJKP)/=ZERO) F_AT(IJKP)=F_AT(IJK)
2618
2619 ENDIF
2620
2621 ENDDO
2622
2623
2624
2625 call send_recv(F_AT,2)
2626
2627
2628
2629 = 0
2630 DO IJK = IJKSTART3, IJKEND3
2631 IF(INTERIOR_CELL_AT(IJK).AND.F_AT(IJK)==UNDEFINED) N_UNDEFINED = N_UNDEFINED + 1
2632 ENDDO
2633
2634 call global_all_sum( N_UNDEFINED, NTOTAL_UNDEFINED )
2635 IF(NTOTAL_UNDEFINED==0) EXIT
2636
2637 ENDDO
2638
2639
2640 call send_recv(F_AT,2)
2641
2642
2643
2644 IF(N_UNDEFINED>0) THEN
2645 WRITE(*,*)'WARNING: UNABLE TO PROPAGATE F_AT ARRAY FROM myPE=.', MyPE
2646 WRITE(*,*)' THIS USUALLY INDICATE A RANK WITH NO ACTIVE CELL'
2647 WRITE(*,*)' YOU MAY NEED TO ADJUST THE GRID PARTITIONNING'
2648 WRITE(*,*)' TO GET BETTER LOAD_BALANCE.'
2649 ENDIF
2650
2651
2652
2653 CASE ('IJK')
2654
2655 DO I=ISTART3,IEND3
2656 DO K=KSTART3,KEND3
2657 F_FOUND = .FALSE.
2658 DO J=JSTART3,JEND3
2659 IJK=FUNIJK(I,J,K)
2660 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2661 F_FOUND = .TRUE.
2662 CURRENT_F = F_AT(IJK)
2663 EXIT
2664 ENDIF
2665 ENDDO
2666 IF(F_FOUND) THEN
2667 DO J=JSTART3,JEND3
2668 IJK=FUNIJK(I,J,K)
2669 IF(F_AT(IJK)==UNDEFINED) THEN
2670 F_AT(IJK)=CURRENT_F
2671 ELSEIF(F_AT(IJK)/=ZERO) THEN
2672 CURRENT_F = F_AT(IJK)
2673 ENDIF
2674 ENDDO
2675 ENDIF
2676 ENDDO
2677 ENDDO
2678
2679 call send_recv(F_AT,2)
2680
2681 DO J=JSTART3,JEND3
2682 DO I=ISTART3,IEND3
2683 F_FOUND = .FALSE.
2684 DO K=KSTART3,KEND3
2685 IJK=FUNIJK(I,J,K)
2686 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2687 F_FOUND = .TRUE.
2688 CURRENT_F = F_AT(IJK)
2689 EXIT
2690 ENDIF
2691 ENDDO
2692 IF(F_FOUND) THEN
2693 DO K=KSTART3,KEND3
2694 IJK=FUNIJK(I,J,K)
2695 IF(F_AT(IJK)==UNDEFINED) THEN
2696 F_AT(IJK)=CURRENT_F
2697 ELSEIF(F_AT(IJK)/=ZERO) THEN
2698 CURRENT_F = F_AT(IJK)
2699 ENDIF
2700 ENDDO
2701 ENDIF
2702 ENDDO
2703 ENDDO
2704
2705 call send_recv(F_AT,2)
2706
2707 DO K=KSTART3,KEND3
2708 DO J=JSTART3,JEND3
2709 F_FOUND = .FALSE.
2710 DO I=ISTART3,IEND3
2711 IJK=FUNIJK(I,J,K)
2712 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2713 F_FOUND = .TRUE.
2714 CURRENT_F = F_AT(IJK)
2715 EXIT
2716 ENDIF
2717 ENDDO
2718 IF(F_FOUND) THEN
2719 DO I=ISTART3,IEND3
2720 IJK=FUNIJK(I,J,K)
2721 IF(F_AT(IJK)==UNDEFINED) THEN
2722 F_AT(IJK)=CURRENT_F
2723 ELSEIF(F_AT(IJK)/=ZERO) THEN
2724 CURRENT_F = F_AT(IJK)
2725 ENDIF
2726 ENDDO
2727 ENDIF
2728 ENDDO
2729 ENDDO
2730
2731 call send_recv(F_AT,2)
2732
2733 CASE ('JKI')
2734
2735 DO J=JSTART3,JEND3
2736 DO I=ISTART3,IEND3
2737 F_FOUND = .FALSE.
2738 DO K=KSTART3,KEND3
2739 IJK=FUNIJK(I,J,K)
2740 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2741 F_FOUND = .TRUE.
2742 CURRENT_F = F_AT(IJK)
2743 EXIT
2744 ENDIF
2745 ENDDO
2746 IF(F_FOUND) THEN
2747 DO K=KSTART3,KEND3
2748 IJK=FUNIJK(I,J,K)
2749 IF(F_AT(IJK)==UNDEFINED) THEN
2750 F_AT(IJK)=CURRENT_F
2751 ELSEIF(F_AT(IJK)/=ZERO) THEN
2752 CURRENT_F = F_AT(IJK)
2753 ENDIF
2754 ENDDO
2755 ENDIF
2756 ENDDO
2757 ENDDO
2758
2759 call send_recv(F_AT,2)
2760
2761 DO K=KSTART3,KEND3
2762 DO J=JSTART3,JEND3
2763 F_FOUND = .FALSE.
2764 DO I=ISTART3,IEND3
2765 IJK=FUNIJK(I,J,K)
2766 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2767 F_FOUND = .TRUE.
2768 CURRENT_F = F_AT(IJK)
2769 EXIT
2770 ENDIF
2771 ENDDO
2772 IF(F_FOUND) THEN
2773 DO I=ISTART3,IEND3
2774 IJK=FUNIJK(I,J,K)
2775 IF(F_AT(IJK)==UNDEFINED) THEN
2776 F_AT(IJK)=CURRENT_F
2777 ELSEIF(F_AT(IJK)/=ZERO) THEN
2778 CURRENT_F = F_AT(IJK)
2779 ENDIF
2780 ENDDO
2781 ENDIF
2782 ENDDO
2783 ENDDO
2784
2785 call send_recv(F_AT,2)
2786
2787 DO I=ISTART3,IEND3
2788 DO K=KSTART3,KEND3
2789 F_FOUND = .FALSE.
2790 DO J=JSTART3,JEND3
2791 IJK=FUNIJK(I,J,K)
2792 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2793 F_FOUND = .TRUE.
2794 CURRENT_F = F_AT(IJK)
2795 EXIT
2796 ENDIF
2797 ENDDO
2798 IF(F_FOUND) THEN
2799 DO J=JSTART3,JEND3
2800 IJK=FUNIJK(I,J,K)
2801 IF(F_AT(IJK)==UNDEFINED) THEN
2802 F_AT(IJK)=CURRENT_F
2803 ELSEIF(F_AT(IJK)/=ZERO) THEN
2804 CURRENT_F = F_AT(IJK)
2805 ENDIF
2806 ENDDO
2807 ENDIF
2808 ENDDO
2809 ENDDO
2810
2811 call send_recv(F_AT,2)
2812
2813
2814 CASE ('KIJ')
2815
2816
2817
2818 DO K=KSTART3,KEND3
2819 DO J=JSTART3,JEND3
2820 F_FOUND = .FALSE.
2821 DO I=ISTART3,IEND3
2822 IJK=FUNIJK(I,J,K)
2823 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2824 F_FOUND = .TRUE.
2825 CURRENT_F = F_AT(IJK)
2826 EXIT
2827 ENDIF
2828 ENDDO
2829 IF(F_FOUND) THEN
2830 DO I=ISTART3,IEND3
2831 IJK=FUNIJK(I,J,K)
2832 IF(F_AT(IJK)==UNDEFINED) THEN
2833 F_AT(IJK)=CURRENT_F
2834 ELSEIF(F_AT(IJK)/=ZERO) THEN
2835 CURRENT_F = F_AT(IJK)
2836 ENDIF
2837 ENDDO
2838 ENDIF
2839 ENDDO
2840 ENDDO
2841
2842 call send_recv(F_AT,2)
2843
2844 DO I=ISTART3,IEND3
2845 DO K=KSTART3,KEND3
2846 F_FOUND = .FALSE.
2847 DO J=JSTART3,JEND3
2848 IJK=FUNIJK(I,J,K)
2849 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2850 F_FOUND = .TRUE.
2851 CURRENT_F = F_AT(IJK)
2852 EXIT
2853 ENDIF
2854 ENDDO
2855 IF(F_FOUND) THEN
2856 DO J=JSTART3,JEND3
2857 IJK=FUNIJK(I,J,K)
2858 IF(F_AT(IJK)==UNDEFINED) THEN
2859 F_AT(IJK)=CURRENT_F
2860 ELSEIF(F_AT(IJK)/=ZERO) THEN
2861 CURRENT_F = F_AT(IJK)
2862 ENDIF
2863 ENDDO
2864 ENDIF
2865 ENDDO
2866 ENDDO
2867
2868
2869 call send_recv(F_AT,2)
2870
2871 DO J=JSTART3,JEND3
2872 DO I=ISTART3,IEND3
2873 F_FOUND = .FALSE.
2874 DO K=KSTART3,KEND3
2875 IJK=FUNIJK(I,J,K)
2876 IF(F_AT(IJK)/=UNDEFINED.AND.F_AT(IJK)/=ZERO) THEN
2877 F_FOUND = .TRUE.
2878 CURRENT_F = F_AT(IJK)
2879 EXIT
2880 ENDIF
2881 ENDDO
2882 IF(F_FOUND) THEN
2883 DO K=KSTART3,KEND3
2884 IJK=FUNIJK(I,J,K)
2885 IF(F_AT(IJK)==UNDEFINED) THEN
2886 F_AT(IJK)=CURRENT_F
2887 ELSEIF(F_AT(IJK)/=ZERO) THEN
2888 CURRENT_F = F_AT(IJK)
2889 ENDIF
2890 ENDDO
2891 ENDIF
2892 ENDDO
2893 ENDDO
2894
2895
2896 call send_recv(F_AT,2)
2897
2898 CASE DEFAULT
2899 IF(myPE == PE_IO) THEN
2900 WRITE(*,*)'CAD_INTERSECT.'
2901 WRITE(*,*)'UNKNOWN CAD_PROPAGATE_ORDER:',CAD_PROPAGATE_ORDER
2902 WRITE(*,*)'ACCEPTABLE VALUES:'
2903 WRITE(*,*)'IJK'
2904 WRITE(*,*)'JKI'
2905 WRITE(*,*)'KIJ'
2906 ENDIF
2907
2908
2909 END SELECT
2910
2911 call send_recv(F_AT,2)
2912
2913 RETURN
2914
2915 END SUBROUTINE CAD_INTERSECT
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929
2930 SUBROUTINE ADD_FACET_AND_SET_BC_ID(IJK,NN)
2931
2932 USE param
2933 USE param1
2934 USE parallel
2935 USE constant
2936 USE run
2937 USE toleranc
2938 USE geometry
2939 USE indices
2940 USE compar
2941 USE sendrecv
2942 USE quadric
2943 USE cutcell
2944 USE polygon
2945 USE stl
2946 USE mpi_utility
2947
2948 IMPLICIT NONE
2949 INTEGER :: IJK,NN
2950
2951 BC_ID(IJK) = BC_ID_STL_FACE(NN)
2952
2953 IF(N_FACET_AT(IJK)<DIM_FACETS_PER_CELL) THEN
2954
2955 N_FACET_AT(IJK) = N_FACET_AT(IJK) + 1
2956 LIST_FACET_AT(IJK,N_FACET_AT(IJK)) = NN
2957
2958 ELSE
2959
2960 WRITE(*,*) ' FATAL ERROR: TOO MANY FACETS IN CELL: ', IJK
2961 CALL MFIX_EXIT(myPE)
2962
2963 ENDIF
2964
2965 END SUBROUTINE ADD_FACET_AND_SET_BC_ID
2966 END MODULE CUT_CELL_PREPROC
2967