MFIX  2016-1
cut_cell_preprocessing.f
Go to the documentation of this file.
1 MODULE cut_cell_preproc
2  USE exit, ONLY: mfix_exit
3  CONTAINS
4 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
5 ! C
6 ! Module name: CUT_CELL_PREPROCESSING C
7 ! Purpose: Perform the cut-cell preprocessing stage: C
8 ! Identify cut cells, define face areas, and volumes C
9 ! Set flags C
10 ! Compute Interpolations factors C
11 ! Compute Non-orthogonality Corrections terms C
12 ! C
13 ! Author: Jeff Dietiker Date: 21-Feb-08 C
14 ! Reviewer: Date: C
15 ! C
16 ! Revision Number # Date: ##-###-## C
17 ! Author: # C
18 ! Purpose: # C
19 ! C
20 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! Gather Data and writes surface(s) defined by all cut cells
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
148 ! C
149 ! Module name: CG_FLOW_TO_VEL C
150 ! Purpose: Convert flow to velocity bc's C
151 ! C
152 ! C
153 ! Author: Jeff Dietiker Date: 21-Feb-08 C
154 ! Reviewer: Date: C
155 ! C
156 ! Revision Number # Date: ##-###-## C
157 ! Author: # C
158 ! Purpose: # C
159 ! C
160 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! G l o b a l P a r a m e t e r s
188 !-----------------------------------------------
189 !-----------------------------------------------
190 ! L o c a l P a r a m e t e r s
191 !-----------------------------------------------
192 !-----------------------------------------------
193 ! L o c a l V a r i a b l e s
194 !-----------------------------------------------
195 !
196 ! loop/variable indices
197  INTEGER :: M, BCV
198 ! Volumetric flow rate computed from mass flow rate
199  DOUBLE PRECISION :: VOLFLOW
200 ! Solids phase volume fraction
201  DOUBLE PRECISION :: EPS
202 ! Average molecular weight
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 ! If gas mass flow is defined convert it to volumetric flow
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 ! If volumetric flow is also specified compare both
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 ! If gas volumetric flow is defined convert it to velocity
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 !Error caught in Check_data_07
259  ENDIF
260  ENDIF
261 
262  ENDIF
263 
264 
265 !
266 ! Do flow conversions for solids phases
267 !
268  DO m = 1, mmax
269 
270  IF(bc_velmag_s(bcv,m)==undefined) THEN
271 !
272 ! If solids mass flow is defined convert it to volumetric flow
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 ! This error will be caught in a previous routine
279  ENDIF
280 !
281 ! If volumetric flow is also specified compare both
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
353 ! C
354 ! Module name: CONVERT_CG_MI_TO_PS C
355 ! Purpose: Convert CG_MI BCs to Point sources C
356 ! C
357 ! C
358 ! Author: Jeff Dietiker Date: 06-Jan-14 C
359 ! Reviewer: Date: C
360 ! C
361 ! Revision Number # Date: ##-###-## C
362 ! Author: # C
363 ! Purpose: # C
364 ! C
365 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! G l o b a l P a r a m e t e r s
394 !-----------------------------------------------
395 !-----------------------------------------------
396 ! L o c a l P a r a m e t e r s
397 !-----------------------------------------------
398 !-----------------------------------------------
399 ! L o c a l V a r i a b l e s
400 !-----------------------------------------------
401 !
402 ! loop/variable indices
403  INTEGER :: IJK, M, NN, BCV
404 !
405  INTEGER :: iproc
406  INTEGER :: NPS,PSV
407 
408 !-----------------------------------------------
409 !
410 
411 ! print*,'Entering test',MyPE
412 #ifdef MPI
413  CALL mpi_barrier(mpi_comm_world,mpierr)
414 #endif
415 
416 ! Each procesor waits for its turn to find cells where to add a point source and updates the list of point sources
417 
418  do iproc = 0,numpes-1
419  if (mype==iproc) Then
420 
421 ! First, find how many point sources are already defined. This could be regular PS from mfix.dat or new ones
422 ! coming from the convertion of CG_MI to PS
423  nps = 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 ! print *,'Last PS=',NPS
431 
432 ! Next loop through all cells, and when a cut-cell with CG_MI is found, add a point source in this cell
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 ! print*,MyPE,NPS
442 
443  ps_defined(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 ! print*,'PS created:',NPS,PS_MASSFLOW_g(NPS),PS_VOLUME(NPS),BC_VOL(BCV)
513  ENDIF
514  ENDIF
515 
516  ENDDO ! IJK Loop
517 
518  endif ! Work done by each processor in same order as rank
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 ! print*,'Leaving test',MyPE
551 ! call mfix_exit(myPE)
552 
553  RETURN
554 
555  END SUBROUTINE convert_cg_mi_to_ps
556 
557 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
558 ! C
559 ! Module name: CONVERT_CG_MI_TO_PS_PE C
560 ! Purpose: Convert CG_MI BCs to Point sources C
561 ! C
562 ! C
563 ! Author: Jeff Dietiker Date: 06-Jan-14 C
564 ! Reviewer: Date: C
565 ! C
566 ! Revision Number # Date: ##-###-## C
567 ! Author: # C
568 ! Purpose: # C
569 ! C
570 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! G l o b a l P a r a m e t e r s
599 !-----------------------------------------------
600 !-----------------------------------------------
601 ! L o c a l P a r a m e t e r s
602 !-----------------------------------------------
603 !-----------------------------------------------
604 ! L o c a l V a r i a b l e s
605 !-----------------------------------------------
606 !
607 ! loop/variable indices
608  INTEGER :: IJK, BCV
609 !
610  INTEGER :: NPS,PSV
611 !
612 !-----------------------------------------------
613 !
614 
615 ! Find the last Point source that is defined. New point sources
616 ! will be added after that.
617 
618 ! print*,'setting bc_type to CG_NSW and exiting'
619 ! DO BCV = 1, DIMENSION_BC
620 ! IF (BC_TYPE_ENUM(BCV) == 'CG_MI') THEN
621 ! BC_TYPE_ENUM(BCV) = 'CG_NSW'
622 ! print*,'Converted CG_MI to CG_FSW for BC#',BCV
623 ! ENDIF
624 ! ENDDO
625 ! RETURN
626 
627  nps = 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 ! read(*,*)
636 
637 ! Loop though each cell. When a CG_MI is found convert it to a single point source
638 ! and change the BC_TYPE to Free-slip
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 ! This is a temporary setting for the solids phase and will need to be generalalized
683  ps_massflow_s(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 ! IF(Normal_S(IJK,2)/=ONE) print*,'Not vertical'
696 ! IF(Normal_S(IJK,2)==ONE) print*,' vertical'
697 
698 ! IF(Normal_S(IJK,2)/=ONE) PS_MASSFLOW_g(NPS) = ZERO
699 
700 ! IF(.NOT.CUT_CELL_AT(IJK)) THEN
701 ! print*,'turn off PS :not a scalar cut cell'
702 ! PS_MASSFLOW_g(NPS) = ZERO
703 ! ENDIF
704 ! IF(.NOT.CUT_U_CELL_AT(IJK)) THEN
705 ! print*,'turn off PS :not a u cut cell'
706 ! PS_MASSFLOW_g(NPS) = ZERO
707 ! ENDIF
708 ! IF(.NOT.CUT_V_CELL_AT(IJK)) THEN
709 ! print*,'turn off PS :not a v cut cell'
710 ! PS_MASSFLOW_g(NPS) = ZERO
711 ! ENDIF
712 ! IF(.NOT.CUT_W_CELL_AT(IJK)) THEN
713 ! print*,'turn off PS :not a w cut cell'
714 ! PS_MASSFLOW_g(NPS) = ZERO
715 ! ENDIF
716 
717 
718 ! PS_MASSFLOW_g(NPS) = ZERO
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 ! ENDIF
724 
725 ! PS_DEFINED(NPS) = .FALSE.
726  ENDIF
727  ENDIF
728  ENDDO
729 
730 ! DO BCV = 1, DIMENSION_BC
731 ! IF (BC_TYPE_ENUM(BCV) == 'CG_MI') THEN
732 ! BC_TYPE_ENUM(BCV) = 'CG_NSW'
733 ! print*,'Converted CG_MI to CG_FSW for BC#',BCV
734 ! ENDIF
735 ! ENDDO
736 
737 100 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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
768 ! C
769 ! Module name: PRINT_CG_HEADER C
770 ! Purpose: Display Cartesian-Grid Header on screen C
771 ! C
772 ! Author: Jeff Dietiker Date: 21-Feb-08 C
773 ! Reviewer: Date: C
774 ! C
775 ! Revision Number # Date: ##-###-## C
776 ! Author: # C
777 ! Purpose: # C
778 ! C
779 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! IF(ADJUST_PROC_DOMAIN_SIZE) THEN
813 ! WRITE(*,10)'EACH PROCESSOR DOMAIN SIZE WILL BE ADJUSTED FOR BETTER LOAD BALANCING.'
814 ! ELSE
815 ! WRITE(*,10)'WARNING: PROCESSOR DOMAIN SIZE WILL BE UNIFORMLY DISTRIBUTED.'
816 ! WRITE(*,10)'THIS COULD RESULT IN VERY POOR LOAD BALANCING.'
817 ! ENDIF
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
833 ! C
834 ! Module name: EVAL_F C
835 ! Purpose: Evaluate the function f(x,y,z) defining the boundary C
836 ! C
837 ! Author: Jeff Dietiker Date: 21-FEB-08 C
838 ! Reviewer: Date: C
839 ! C
840 ! C
841 ! Literature/Document References: C
842 ! C
843 ! Variables referenced: C
844 ! Variables modified: C
845 ! C
846 ! Local variables: C
847 ! C
848 ! Modified: ## Date: ##-###-## C
849 ! Purpose: ## C
850 ! C
851 ! Literature/Document References: ## C
852 ! C
853 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! CLIP_FLAG=.FALSE.
900  CALL get_f_quadric(x1,x2,x3,q_id,f_g(group,0),clip_flag)
901 ! CLIP_FLAG=.TRUE.
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 ! CASE('STL')
930 
931 ! CALL EVAL_STL_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
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 ! WRITE(*,*)'STL'
942  CALL mfix_exit(mype)
943 
944  END SELECT
945 
946  DEALLOCATE(f_g)
947 
948  RETURN
949  END SUBROUTINE eval_f
950 
951 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
952 ! C
953 ! Module name: intersect_line C
954 ! Purpose: Finds the intersection between the quadric surface , C
955 ! and the line (xa,ya,za) and (xb,yb,zb). C
956 ! C
957 ! Author: Jeff Dietiker Date: 21-Feb-08 C
958 ! Reviewer: Date: C
959 ! C
960 ! Revision Number # Date: ##-###-## C
961 ! Author: # C
962 ! Purpose: # C
963 ! C
964 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! Initial guesses
990  y1 = 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 ! The line from (x1,y1,z1) and (x2,y2,z2) is parametrized
1004 ! from t1 = ZERO to t2 = ONE
1005 !======================================================================
1006 
1007  niter = 1
1008 
1009  clip_flag = (clip_flag1).AND.(clip_flag2)
1010 
1011  if(dabs(f1)<tol_f) then ! ignore intersection at corner
1012  xc = undefined
1013  yc = undefined
1014  zc = undefined
1015  intersect_flag = .false.
1016  elseif(dabs(f2)<tol_f) then ! ignore intersection at corner
1017  xc = 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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1084 ! C
1085 ! Module name: INTERSECT C
1086 ! Purpose: Intersects quadric with grid C
1087 ! C
1088 ! Author: Jeff Dietiker Date: 21-Feb-08 C
1089 ! Reviewer: Date: C
1090 ! C
1091 ! Revision Number # Date: ##-###-## C
1092 ! Author: # C
1093 ! Purpose: # C
1094 ! C
1095 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! Get coordinates of eight nodes
1126 !======================================================================
1127 
1128  CALL get_cell_node_coordinates(ijk,type_of_cell)
1129 
1130 !======================================================================
1131 ! Intersection with Edge 7 (node 7-8, Face North-Top):
1132 !======================================================================
1133  xa = 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 ! IF(USE_STL) THEN
1170 ! CALL INTERSECT_LINE_WITH_STL(xa,ya,za,xb,yb,zb,INTERSECT_X(IJK),xc,yc,zc)
1171 ! IF(INTERSECT_X(IJK)) Xi = xc
1172 ! ENDIF
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 ! Intersection with Edge 6 (node 6-8, Face East-Top):
1186 !======================================================================
1187  xa = 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 ! IF(USE_STL) THEN
1218 ! CALL INTERSECT_LINE_WITH_STL(xa,ya,za,xb,yb,zb,INTERSECT_Y(IJK),xc,yc,zc)
1219 ! IF(INTERSECT_Y(IJK)) Yi = yc
1220 ! ENDIF
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 ! Intersection with Edge 11 (node 4-8, Face East-North):
1233 !======================================================================
1234  xa = 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  ! IF(USE_STL) THEN
1265  ! CALL INTERSECT_LINE_WITH_STL(xa,ya,za,xb,yb,zb,INTERSECT_Z(IJK),xc,yc,zc)
1266  ! IF(INTERSECT_Z(IJK)) Zi = zc
1267  ! ENDIF
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1303 ! C
1304 ! Module name: CLEAN_INTERSECT C
1305 ! Purpose: Remove Intersection flags in preparation of small cell C
1306 ! removal C
1307 ! C
1308 ! Author: Jeff Dietiker Date: 21-Feb-08 C
1309 ! Reviewer: Date: C
1310 ! C
1311 ! Revision Number # Date: ##-###-## C
1312 ! Author: # C
1313 ! Purpose: # C
1314 ! C
1315 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! When inputing geometry from CAD (STL or MSH file), the snapping procedure is
1346 ! dependent on the value of F at the cell corners
1347 ! For other gemoetry inputs (say quadrics), This is not needed, and the value
1348 ! of F_TEST is set to .TRUE. here
1349  cad = use_msh.OR.use_stl
1350  f_test = .true.
1351 
1352 !======================================================================
1353 ! Get coordinates of eight nodes
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 ! Clean Intersection with Edge 7 (node 7-8, Face North-Top):
1405 !======================================================================
1406 
1407  xa = 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) ! MAXIMUM DISTANCE FROM CORNER
1416 
1417  IF(intersect_x(ijk)) THEN
1418 
1419  dfc = dabs(xi-xa) ! DISTANCE FROM CORNER (NODE 7)
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) ! DISTANCE FROM CORNER (NODE 8)
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 ! Clean Intersection with Edge 6 (node 6-8, Face East-Top):
1477 !======================================================================
1478 
1479  xa = 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) ! MAXIMUM DISTANCE FROM CORNER
1484 
1485 
1486  IF(intersect_y(ijk)) THEN
1487 
1488  dfc = dabs(yi-ya) ! DISTANCE FROM CORNER (NODE 6)
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) ! DISTANCE FROM CORNER (NODE 8)
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 ! Intersection with Edge 11 (node 4-8, Face East-North):
1545 !======================================================================
1546 
1547  xa = 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) ! MAXIMUM DISTANCE FROM CORNER
1552 
1553  IF(intersect_z(ijk)) THEN
1554 
1555  dfc = dabs(zi-za) ! DISTANCE FROM CORNER (NODE 4)
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) ! DISTANCE FROM CORNER (NODE 8)
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 ! F_AT(IJKM) = UNDEFINED
1593 ! IF(F_AT(IJK)/=ZERO) SNAP(IJK) = .TRUE.
1594 ! F_AT(IJK) = ZERO
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1620 ! C
1621 ! Module name: CLEAN_INTERSECT_SCALAR C
1622 ! Purpose: Remove Intersection flags in preparation of small cell C
1623 ! removal C
1624 ! C
1625 ! Author: Jeff Dietiker Date: 04-Dec-14 C
1626 ! Reviewer: Date: C
1627 ! C
1628 ! Revision Number # Date: ##-###-## C
1629 ! Author: # C
1630 ! Purpose: # C
1631 ! C
1632 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1692 ! C
1693 ! Module name: SET_SNAP_FLAG C
1694 ! Purpose: Set SNAP flag in preparation of intersection flag removal C
1695 ! removal C
1696 ! C
1697 ! Author: Jeff Dietiker Date: 04-Dec-14 C
1698 ! Reviewer: Date: C
1699 ! C
1700 ! Revision Number # Date: ##-###-## C
1701 ! Author: # C
1702 ! Purpose: # C
1703 ! C
1704 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! When inputing geometry from CAD (STL or MSH file), the snapping procedure is
1735 ! dependent on the value of F at the cell corners
1736 ! For other gemoetry inputs (say quadrics), This is not needed, and the value
1737 ! of F_TEST is set to .TRUE. here
1738  cad = use_msh.OR.use_stl
1739  f_test = .true.
1740 
1741 !======================================================================
1742 ! Get coordinates of eight nodes
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 ! Clean Intersection with Edge 7 (node 7-8, Face North-Top):
1765 !======================================================================
1766 
1767  xa = 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) ! MAXIMUM DISTANCE FROM CORNER
1777 
1778  IF(intersect_x(ijk)) THEN
1779 
1780  dfc = dabs(xi-xa) ! DISTANCE FROM CORNER (NODE 7)
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) ! DISTANCE FROM CORNER (NODE 8)
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 ! Clean Intersection with Edge 6 (node 6-8, Face East-Top):
1818 !======================================================================
1819 
1820  xa = 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) ! MAXIMUM DISTANCE FROM CORNER
1825 
1826 
1827  IF(intersect_y(ijk)) THEN
1828 
1829  dfc = dabs(yi-ya) ! DISTANCE FROM CORNER (NODE 6)
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) ! DISTANCE FROM CORNER (NODE 8)
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 ! Intersection with Edge 11 (node 4-8, Face East-North):
1870 !======================================================================
1871 
1872  xa = 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) ! MAXIMUM DISTANCE FROM CORNER
1877 
1878  IF(intersect_z(ijk)) THEN
1879 
1880  dfc = dabs(zi-za) ! DISTANCE FROM CORNER (NODE 4)
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) ! DISTANCE FROM CORNER (NODE 8)
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  !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1923 ! C
1924 ! Module name: REMOVE_INTERSECT_FLAG C
1925 ! Purpose: Remove Intersection flags C
1926 ! removal C
1927 ! C
1928 ! Author: Jeff Dietiker Date: 04-Dec-14 C
1929 ! Reviewer: Date: C
1930 ! C
1931 ! Revision Number # Date: ##-###-## C
1932 ! Author: # C
1933 ! Purpose: # C
1934 ! C
1935 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1990 ! C
1991 ! Module name: OPEN_CUT_CELL_FILES C
1992 ! Purpose: Open CUT CELL related file C
1993 ! and writes headers C
1994 ! C
1995 ! Author: Jeff Dietiker Date: 21-Feb-08 C
1996 ! Reviewer: Date: C
1997 ! C
1998 ! Revision Number # Date: ##-###-## C
1999 ! Author: # C
2000 ! Purpose: # C
2001 ! C
2002 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2020 ! C
2021 ! Module name: CLOSE_CUT_CELL_FILES C
2022 ! Purpose: Close CUT CELL related file C
2023 ! C
2024 ! Author: Jeff Dietiker Date: 21-Feb-08 C
2025 ! Reviewer: Date: C
2026 ! C
2027 ! Revision Number # Date: ##-###-## C
2028 ! Author: # C
2029 ! Purpose: # C
2030 ! C
2031 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2051 ! C
2052 ! Module name: CAD_INTERSECT C
2053 ! Purpose: Intersects CAD (STL file or MSH file) geometry with grid C
2054 ! C
2055 ! Author: Jeff Dietiker Date: 21-Feb-08 C
2056 ! Reviewer: Date: C
2057 ! C
2058 ! Revision Number # Date: ##-###-## C
2059 ! Author: # C
2060 ! Purpose: # C
2061 ! C
2062 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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 ! CHARACTER (LEN=3) :: CAD_PROPAGATE_ORDER
2107 
2108  intersect_x = .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 ! Get coordinates of eight nodes
2266 !======================================================================
2267 
2268  CALL get_cell_node_coordinates(ijk,type_of_cell)
2269 
2270 !======================================================================
2271 ! Intersection with Edge 7 (node 7-8, Face North-Top):
2272 !======================================================================
2273  xa = 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 ! Check if intersection occurs at corners
2282 
2283  CALL is_point_inside_facet(xa,ya,za,nn,inside_facet_a)
2284 
2285  IF(inside_facet_a) THEN ! corner intersection at node 7
2286 
2287  f_at(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 ! corner intersection at node 8
2297 
2298  f_at(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 ! Check intersection within line 7-8, excluding corners
2307 
2308  intersect_flag = .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  ! Ignore intersections when two intersections are detected on the same edge
2319  intersect_x(ijk) = .false.
2320  ENDIF
2321  ELSE
2322  intersect_x(ijk) = .true.
2323  xint(ijk) = xc
2324 
2325 ! Set values at corners if they are not zero
2326 
2327  n7(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 ! Intersection with Edge 6 (node 6-8, Face East-Top):
2356 !======================================================================
2357  xa = x_node(6)
2358  ya = y_node(6)
2359  za = z_node(6)
2360 
2361 ! Check if intersection occurs at corners
2362 
2363  CALL is_point_inside_facet(xa,ya,za,nn,inside_facet_a)
2364 
2365  IF(inside_facet_a) THEN ! corner intersection at node 6
2366 
2367  f_at(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 ! Check intersection within line 6-8, excluding corners
2376 
2377 
2378 
2379  intersect_flag = .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. ! Ignore intersections when two intersections are detected on the same edge
2393 
2394  ENDIF
2395 
2396  ELSE
2397 
2398 
2399  intersect_y(ijk) = .true.
2400  yint(ijk) = yc
2401 
2402 ! Set values at corners if they are not zero
2403 
2404  n6(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 ! Intersection with Edge 11 (node 4-8, Face East-North):
2436 !======================================================================
2437  xa = x_node(4)
2438  ya = y_node(4)
2439  za = z_node(4)
2440 
2441 ! Check if intersection occurs at corners
2442 
2443  CALL is_point_inside_facet(xa,ya,za,nn,inside_facet_a)
2444 
2445  IF(inside_facet_a) THEN ! corner intersection at node 4
2446 
2447  f_at(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 ! Check intersection within line 4-8, excluding corners
2457 
2458 
2459  intersect_flag = .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. ! Ignore intersections when two intersections are detected on the same edge
2472 
2473  ENDIF
2474 
2475  ELSE
2476 
2477 
2478  intersect_z(ijk) = .true.
2479  zint(ijk) = zc
2480 
2481 
2482 ! Set values at corners if they are not zero
2483 
2484  n4(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 ! I loop
2519  ENDDO ! J loop
2520  ENDDO ! K loop
2521 
2522 
2523  ENDDO ! Loop over facets
2524 
2525  current_f = undefined
2526 
2527 ! Overwrite small values to set them to zero
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 ! Propagates node values to all interior cells
2536 ! in order defined by CAD_PROPAGATE_ORDER (outer loop)
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 ! Clean-up intersection flags in preparaton of small cells removal
2548 !======================================================================
2549 
2550  DO ijk = ijkstart3, ijkend3
2551 
2552 ! IF(INTERIOR_CELL_AT(IJK)) THEN
2553 
2554 ! IF(POTENTIAL_CUT_CELL_AT(IJK)) CALL CLEAN_INTERSECT(IJK,'SCALAR',Xn_int(IJK),Ye_int(IJK),Zt_int(IJK))
2555 
2556  CALL clean_intersect(ijk,type_of_cell,xint(ijk),yint(ijk),zint(ijk))
2557 
2558 ! ENDIF
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 ! After intersecting the edges of the background mesh with the STL facets,
2573 ! the end points (i.e., cell corners) are assigned a value, called F_AT, where:
2574 ! F_AT = zero if the corner point is on a facet (within some tolerance TOL_STL),
2575 ! F_AT < zero if the corner point is inside the fluid region,
2576 ! F_AT > zero if the corner point is outside the fluid region.
2577 ! At this point F_AT is only defined across edges that intersect the STL facets,
2578 ! and it must be propagated to all cell corners to determine if uncut cells
2579 ! are inside or outside the fluid region.
2580 
2581 ! Only F_AT values that are defined and not zero are propagated to their direct
2582 ! neighbors, if it is not already defined. The propagation is repeated
2583 ! at most N_PROPMAX. The loop is exited when all F_AT values are defined.
2584 ! N_PROPMAX could be increased for very large domains.
2585 ! The propagation of F_AT will stop anytime a boundary is encountered since F_AT
2586 ! changes sign across a boundary.
2587 !
2588  DO n_prop=1,n_propmax
2589 
2590  DO ijk = ijkstart3, ijkend3
2591 
2592 ! Aaron
2593  i = 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 ! End aaron
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 ! IJK Loop
2622 
2623 
2624 ! Communicate F_AT accross processors for DMP runs
2625  call send_recv(f_at,2)
2626 
2627 ! Count the number of undefined values of F_AT
2628 ! and exit loop if all values of F_AT have been propagated
2629  n_undefined = 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 ! N_PROP Loop
2638 
2639 
2640  call send_recv(f_at,2)
2641 
2642 ! If a process still has undefined values of F_AT, this probably means
2643 ! that all cells belonging to that process are dead cells.
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 ! CALL MFIX_EXIT(myPE)
2908 
2909  END SELECT
2910 
2911  call send_recv(f_at,2)
2912 
2913  RETURN
2914 
2915  END SUBROUTINE cad_intersect
2916 
2917 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2918 ! C
2919 ! Module name: ADD_FACET C
2920 ! Purpose: Add facet to list in IJK scalar cell C
2921 ! C
2922 ! Author: Jeff Dietiker Date: 15-Oct-13 C
2923 ! Reviewer: Date: C
2924 ! C
2925 ! Revision Number # Date: ##-###-## C
2926 ! Author: # C
2927 ! Purpose: # C
2928 ! C
2929 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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) ! Set tentative BC_ID
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
stl::n_facets
integer n_facets
Definition: stl_mod.f:8
ps::ps_i_w
integer, dimension(dimension_ps) ps_i_w
Definition: ps_mod.f:40
mpi_utility
Definition: mpi_utility_mod.f:6
compar::iend3
integer iend3
Definition: compar_mod.f:80
cutcell::re_indexing
logical re_indexing
Definition: cutcell_mod.f:16
ps::ps_v_g
double precision, dimension(dimension_ps) ps_v_g
Definition: ps_mod.f:52
funits::dmp_log
logical dmp_log
Definition: funits_mod.f:6
set_3d_cut_w_cell_flags
subroutine set_3d_cut_w_cell_flags
Definition: get_cut_cell_flags.f:1112
bc::bc_volflow_g
double precision, dimension(dimension_bc) bc_volflow_g
Definition: bc_mod.f:195
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
cutcell::potential_cut_cell_at
logical, dimension(:), allocatable potential_cut_cell_at
Definition: cutcell_mod.f:474
ps::ps_t_g
double precision, dimension(dimension_ps) ps_t_g
Definition: ps_mod.f:62
ps::ps_v_s
double precision, dimension(dimension_ps, dim_m) ps_v_s
Definition: ps_mod.f:70
compar::jstart3
integer jstart3
Definition: compar_mod.f:80
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
eval_stl_fct_at
subroutine eval_stl_fct_at(TYPE_OF_CELL, IJK, NODE, f_stl, CLIP_FLAG, B
Definition: get_stl_data.f:1085
cutcell::yg_n
double precision, dimension(:), allocatable yg_n
Definition: cutcell_mod.f:45
cutcell::xn_int
double precision, dimension(:), allocatable xn_int
Definition: cutcell_mod.f:333
bc::bc_t_g
double precision, dimension(dimension_bc) bc_t_g
Definition: bc_mod.f:97
cdist::bdist_io
logical bdist_io
Definition: cdist_mod.f:4
stl::list_facet_at
integer, dimension(:,:), allocatable list_facet_at
Definition: stl_mod.f:57
toleranc::compare
logical function compare(V1, V2)
Definition: toleranc_mod.f:94
set_3d_cut_u_cell_flags
subroutine set_3d_cut_u_cell_flags
Definition: get_cut_cell_flags.f:653
cut_cell_preproc::set_snap_flag
subroutine set_snap_flag(IJK, TYPE_OF_CELL, Xi, Yi, Zi)
Definition: cut_cell_preprocessing.f:1706
cut_cell_preproc::open_cut_cell_files
subroutine open_cut_cell_files
Definition: cut_cell_preprocessing.f:2004
quadric::itermax_int
integer itermax_int
Definition: quadric_mod.f:78
quadric::group_size
integer, dimension(dim_group) group_size
Definition: quadric_mod.f:68
stl::dim_facets_per_cell
integer dim_facets_per_cell
Definition: stl_mod.f:55
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
get_w_master_cells
subroutine get_w_master_cells
Definition: get_master.f:313
cutcell::xg_e
double precision, dimension(:), allocatable xg_e
Definition: cutcell_mod.f:44
cut_cell_preproc::convert_cg_mi_to_ps
subroutine convert_cg_mi_to_ps
Definition: cut_cell_preprocessing.f:367
vtk::write_vtk_files
logical write_vtk_files
Definition: vtk_mod.f:24
cutcell::ye_int
double precision, dimension(:), allocatable ye_int
Definition: cutcell_mod.f:338
stl::vertex
double precision, dimension(3, 3, dim_stl) vertex
Definition: stl_mod.f:20
cut_cell_preproc::cad_intersect
subroutine cad_intersect(TYPE_OF_CELL, Xint, Yint, Zint)
Definition: cut_cell_preprocessing.f:2064
cutcell::z_node
double precision, dimension(0:15) z_node
Definition: cutcell_mod.f:76
cut_cell_preproc::cg_flow_to_vel
subroutine cg_flow_to_vel
Definition: cut_cell_preprocessing.f:162
is_point_inside_facet
subroutine is_point_inside_facet(Px, Py, Pz, FACET, INSIDE_FACET)
Definition: get_stl_data.f:1320
gather_data
subroutine gather_data
Definition: vtk_out.f:1700
bc::bc_w_s
double precision, dimension(dimension_bc, dim_m) bc_w_s
Definition: bc_mod.f:129
get_3d_alpha_u_cut_cell
subroutine get_3d_alpha_u_cut_cell
Definition: get_alpha.f:16
send_receive_1d_logical
subroutine send_receive_1d_logical(L1D, NLAYERS)
Definition: dmp_cartesian.f:277
cutcell::cg_safe_mode
integer, dimension(10) cg_safe_mode
Definition: cutcell_mod.f:415
cutcell::y_node
double precision, dimension(0:15) y_node
Definition: cutcell_mod.f:75
bc::bc_x_s
double precision, dimension(dimension_bc, dim_m, dim_n_s) bc_x_s
Definition: bc_mod.f:254
get_3d_alpha_v_cut_cell
subroutine get_3d_alpha_v_cut_cell
Definition: get_alpha.f:421
set_3d_cut_cell_treatment_flags
subroutine set_3d_cut_cell_treatment_flags
Definition: get_cut_cell_flags.f:1318
cutcell::print_warnings
logical print_warnings
Definition: cutcell_mod.f:416
get_f_quadric
subroutine get_f_quadric(x1, x2, x3, Q_ID, f, CLIP_FLAG)
Definition: define_quadrics.f:83
stl::n_facet_at
integer, dimension(:), allocatable n_facet_at
Definition: stl_mod.f:56
bc::bc_vol
double precision, dimension(dimension_bc) bc_vol
Definition: bc_mod.f:248
get_3d_alpha_w_cut_cell
subroutine get_3d_alpha_w_cut_cell
Definition: get_alpha.f:827
set_odxyz_v_cut_cell
subroutine set_odxyz_v_cut_cell
Definition: set_Odxyz.f:147
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
compar::kstart3
integer kstart3
Definition: compar_mod.f:80
ps::ps_j_n
integer, dimension(dimension_ps) ps_j_n
Definition: ps_mod.f:43
param::dimension_bc
integer, parameter dimension_bc
Definition: param_mod.f:61
vtk
Definition: vtk_mod.f:1
stl::cad_propagate_order
character(len=3) cad_propagate_order
Definition: stl_mod.f:77
get_u_master_cells
subroutine get_u_master_cells
Definition: get_master.f:15
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
cutcell::normal_s
double precision, dimension(:,:), allocatable normal_s
Definition: cutcell_mod.f:120
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
compar::mpierr
integer mpierr
Definition: compar_mod.f:27
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
set_3d_cut_cell_flags
subroutine set_3d_cut_cell_flags
Definition: get_cut_cell_flags.f:16
physprop::calc_mw
double precision function calc_mw(X_G, DIM, L, NMAX, MW_G)
Definition: physprop_mod.f:176
quadric::relation_with_previous
character(len=9), dimension(dim_group) relation_with_previous
Definition: quadric_mod.f:74
ps::ps_massflow_g
double precision, dimension(dimension_ps) ps_massflow_g
Definition: ps_mod.f:48
cutcell::cg_header_was_printed
logical cg_header_was_printed
Definition: cutcell_mod.f:465
ps::ps_defined
logical, dimension(dimension_ps) ps_defined
Definition: ps_mod.f:29
set_3d_cut_v_cell_flags
subroutine set_3d_cut_v_cell_flags
Definition: get_cut_cell_flags.f:884
stl::bc_id_stl_face
integer, dimension(dim_stl) bc_id_stl_face
Definition: stl_mod.f:51
physprop::mw_g
double precision, dimension(dim_n_g) mw_g
Definition: physprop_mod.f:124
eval_usr_fct
subroutine eval_usr_fct(x1, x2, x3, Q, f_usr, CLIP_FLAG)
Definition: eval_usr_fct.f:20
cut_cell_preproc::add_facet_and_set_bc_id
subroutine add_facet_and_set_bc_id(IJK, NN)
Definition: cut_cell_preprocessing.f:2931
ps::ps_w_g
double precision, dimension(dimension_ps) ps_w_g
Definition: ps_mod.f:53
bc::bc_v_g
double precision, dimension(dimension_bc) bc_v_g
Definition: bc_mod.f:117
bc::bc_velmag_s
double precision, dimension(dimension_bc, dim_m) bc_velmag_s
Definition: bc_mod.f:140
ps::ps_t_s
double precision, dimension(dimension_ps, dim_m) ps_t_s
Definition: ps_mod.f:80
compar::kend3
integer kend3
Definition: compar_mod.f:80
compar::numpes
integer numpes
Definition: compar_mod.f:24
quadric::n_group
integer n_group
Definition: quadric_mod.f:66
ps::ps_volume
double precision, dimension(dimension_ps) ps_volume
Definition: ps_mod.f:84
cut_cell_preproc::clean_intersect_scalar
subroutine clean_intersect_scalar
Definition: cut_cell_preprocessing.f:1634
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
compar::pe_io
integer pe_io
Definition: compar_mod.f:30
bc::bc_volflow_s
double precision, dimension(dimension_bc, dim_m) bc_volflow_s
Definition: bc_mod.f:198
cutcell::f_at
double precision, dimension(:), allocatable f_at
Definition: cutcell_mod.f:486
ps::ps_k_b
integer, dimension(dimension_ps) ps_k_b
Definition: ps_mod.f:44
ps::ps_u_s
double precision, dimension(dimension_ps, dim_m) ps_u_s
Definition: ps_mod.f:69
cut_cell_preproc::remove_intersect_flag
subroutine remove_intersect_flag(IJK)
Definition: cut_cell_preprocessing.f:1937
sendrecv::send_recv
Definition: sendrecv_mod.f:75
geometry::kmax1
integer kmax1
Definition: geometry_mod.f:58
eos::eosg
double precision function eosg(MW, PG, TG)
Definition: eos_mod.f:22
cg_get_bc_area
subroutine cg_get_bc_area
Definition: get_bc_area.f:129
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
cutcell::use_stl
logical use_stl
Definition: cutcell_mod.f:428
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
bc::bc_t_s
double precision, dimension(dimension_bc, dim_m) bc_t_s
Definition: bc_mod.f:101
physprop::ro_g0
double precision ro_g0
Definition: physprop_mod.f:59
cut_cell_preproc::intersect_line
subroutine intersect_line(METHOD, xa, ya, za, xb, yb, zb, Q_ID, INTERSECT_FLAG, xc, yc, zc)
Definition: cut_cell_preprocessing.f:966
geometry::imax1
integer imax1
Definition: geometry_mod.f:54
compar::jend3
integer jend3
Definition: compar_mod.f:80
stl
Definition: stl_mod.f:1
exit
Definition: exit.f:2
cdist
Definition: cdist_mod.f:2
get_cell_node_coordinates
subroutine get_cell_node_coordinates(IJK, TYPE_OF_CELL)
Definition: get_connectivity.f:324
eos
Definition: eos_mod.f:10
set_odxyz_w_cut_cell
subroutine set_odxyz_w_cut_cell
Definition: set_Odxyz.f:282
print_grid_statistics
subroutine print_grid_statistics
Definition: vtk_out.f:2129
cut_cell_preproc::print_cg_header
subroutine print_cg_header
Definition: cut_cell_preprocessing.f:781
stl::tol_stl
double precision tol_stl
Definition: stl_mod.f:45
stl::norm_face
double precision, dimension(3, dim_stl) norm_face
Definition: stl_mod.f:22
quadric::tol_f
double precision tol_f
Definition: quadric_mod.f:76
bc::bc_p_g
double precision, dimension(dimension_bc) bc_p_g
Definition: bc_mod.f:80
cutcell::zt_int
double precision, dimension(:), allocatable zt_int
Definition: cutcell_mod.f:343
geometry::dx
double precision, dimension(0:dim_i) dx
Definition: geometry_mod.f:68
cutcell::n_usr_def
integer n_usr_def
Definition: cutcell_mod.f:424
cutcell::tol_snap
double precision, dimension(3) tol_snap
Definition: cutcell_mod.f:374
mpi_utility::global_all_sum
Definition: mpi_utility_mod.f:66
cutcell::intersect_z
logical, dimension(:), allocatable intersect_z
Definition: cutcell_mod.f:66
ps::ps_k_t
integer, dimension(dimension_ps) ps_k_t
Definition: ps_mod.f:45
geometry::xlength
double precision xlength
Definition: geometry_mod.f:33
param1::half
double precision, parameter half
Definition: param1_mod.f:28
geometry::jmax1
integer jmax1
Definition: geometry_mod.f:56
funits::unit_log
integer, parameter unit_log
Definition: funits_mod.f:21
run
Definition: run_mod.f:13
cut_cell_preproc::cut_cell_preprocessing
subroutine cut_cell_preprocessing
Definition: cut_cell_preprocessing.f:22
bc::bc_velmag_g
double precision, dimension(dimension_bc) bc_velmag_g
Definition: bc_mod.f:136
ps::ps_x_g
double precision, dimension(dimension_ps, dim_n_g) ps_x_g
Definition: ps_mod.f:59
param
Definition: param_mod.f:2
define_quadrics
subroutine define_quadrics
Definition: define_quadrics.f:15
param::dimension_ps
integer, parameter dimension_ps
Definition: param_mod.f:65
cut_cell_preproc::close_cut_cell_files
subroutine close_cut_cell_files
Definition: cut_cell_preprocessing.f:2033
bc::bc_v_s
double precision, dimension(dimension_bc, dim_m) bc_v_s
Definition: bc_mod.f:121
flow_to_vel
subroutine flow_to_vel(DO_VEL_CHECK)
Definition: flow_to_vel.f:558
bc::bc_massflow_g
double precision, dimension(dimension_bc) bc_massflow_g
Definition: bc_mod.f:201
cutcell::interior_cell_at
logical, dimension(:), allocatable interior_cell_at
Definition: cutcell_mod.f:40
physprop::nmax
integer, dimension(0:dim_m) nmax
Definition: physprop_mod.f:119
geometry::jmin1
integer jmin1
Definition: geometry_mod.f:42
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
cutcell::cut_cell_at
logical, dimension(:), allocatable cut_cell_at
Definition: cutcell_mod.f:355
bc::bc_massflow_s
double precision, dimension(dimension_bc, dim_m) bc_massflow_s
Definition: bc_mod.f:204
compar::mype
integer mype
Definition: compar_mod.f:24
bc::bc_u_g
double precision, dimension(dimension_bc) bc_u_g
Definition: bc_mod.f:109
cut_cell_preproc::clean_intersect
subroutine clean_intersect(IJK, TYPE_OF_CELL, Xi, Yi, Zi)
Definition: cut_cell_preprocessing.f:1317
physprop::mw_avg
double precision mw_avg
Definition: physprop_mod.f:71
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
cutcell::use_msh
logical use_msh
Definition: cutcell_mod.f:430
set_ghost_cell_flags
subroutine set_ghost_cell_flags
Definition: get_cut_cell_flags.f:1449
ps::ps_massflow_s
double precision, dimension(dimension_ps, dim_m) ps_massflow_s
Definition: ps_mod.f:66
allocate_cut_cell_arrays
subroutine allocate_cut_cell_arrays
Definition: allocate_cut_cell_arrays.f:3
bc::bc_u_s
double precision, dimension(dimension_bc, dim_m) bc_u_s
Definition: bc_mod.f:113
polygon::n_polygon
integer n_polygon
Definition: polygon_mod.f:6
eval_poly_fct
subroutine eval_poly_fct(x1, x2, x3, Q, f_pol, CLIP_FLAG, BCID)
Definition: get_poly_data.f:132
quadric::group_relation
character(len=9), dimension(dim_group) group_relation
Definition: quadric_mod.f:72
scales::p_ref
double precision p_ref
Definition: scales_mod.f:10
cutcell::bc_id
integer, dimension(:), allocatable bc_id
Definition: cutcell_mod.f:433
ps::ps_w_s
double precision, dimension(dimension_ps, dim_m) ps_w_s
Definition: ps_mod.f:71
cutcell::intersect_x
logical, dimension(:), allocatable intersect_x
Definition: cutcell_mod.f:64
bc::bc_x_g
double precision, dimension(dimension_bc, dim_n_g) bc_x_g
Definition: bc_mod.f:251
bc::bc_ep_g
double precision, dimension(dimension_bc) bc_ep_g
Definition: bc_mod.f:77
get_distance_to_wall
subroutine get_distance_to_wall
Definition: get_delh.f:592
ps::ps_x_s
double precision, dimension(dimension_ps, dim_m, dim_n_s) ps_x_s
Definition: ps_mod.f:77
ps
Definition: ps_mod.f:22
cut_cell_preproc::intersect
subroutine intersect(IJK, TYPE_OF_CELL, Xi, Yi, Zi)
Definition: cut_cell_preprocessing.f:1097
physprop::ro_s0
double precision, dimension(dim_m) ro_s0
Definition: physprop_mod.f:28
quadric::dim_quadric
integer, parameter dim_quadric
Definition: quadric_mod.f:4
mpi_utility::bcast
Definition: mpi_utility_mod.f:53
cut_cell_preproc
Definition: cut_cell_preprocessing.f:1
cutcell::intersect_y
logical, dimension(:), allocatable intersect_y
Definition: cutcell_mod.f:65
geometry::ylength
double precision ylength
Definition: geometry_mod.f:35
cut_cell_preproc::eval_f
subroutine eval_f(METHOD, x1, x2, x3, Q, f, CLIP_FLAG)
Definition: cut_cell_preprocessing.f:855
ps::point_source
logical point_source
Definition: ps_mod.f:27
write_cut_surface_vtk
subroutine write_cut_surface_vtk
Definition: vtk_out.f:1470
geometry::imin1
integer imin1
Definition: geometry_mod.f:40
set_odxyz_u_cut_cell
subroutine set_odxyz_u_cut_cell
Definition: set_Odxyz.f:16
send_receive_cut_cell_variables
subroutine send_receive_cut_cell_variables
Definition: dmp_cartesian.f:15
cut_cell_preproc::convert_cg_mi_to_ps_pe
subroutine convert_cg_mi_to_ps_pe
Definition: cut_cell_preprocessing.f:572
ps::ps_u_g
double precision, dimension(dimension_ps) ps_u_g
Definition: ps_mod.f:51
ps::ps_j_s
integer, dimension(dimension_ps) ps_j_s
Definition: ps_mod.f:42
geometry::vol
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
compar::istart3
integer istart3
Definition: compar_mod.f:80
cutcell::zg_t
double precision, dimension(:), allocatable zg_t
Definition: cutcell_mod.f:46
bc::bc_w_g
double precision, dimension(dimension_bc) bc_w_g
Definition: bc_mod.f:125
bc::cg_mi_converted_to_ps
logical, dimension(dimension_bc) cg_mi_converted_to_ps
Definition: bc_mod.f:406
quadric::group_q
integer, dimension(dim_group, dim_quadric) group_q
Definition: quadric_mod.f:70
ps::ps_i_e
integer, dimension(dimension_ps) ps_i_e
Definition: ps_mod.f:41
geometry::kmin1
integer kmin1
Definition: geometry_mod.f:44
geometry::x
double precision, dimension(:), allocatable x
Definition: geometry_mod.f:129
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
geometry::zlength
double precision zlength
Definition: geometry_mod.f:37
cutcell::x_node
double precision, dimension(0:15) x_node
Definition: cutcell_mod.f:74
intersect_line_with_facet
subroutine intersect_line_with_facet(xa, ya, za, xb, yb, zb, FACET, INTERSECT
Definition: get_stl_data.f:1175
bc::bc_rop_s
double precision, dimension(dimension_bc, dim_m) bc_rop_s
Definition: bc_mod.f:92
bc
Definition: bc_mod.f:23
get_v_master_cells
subroutine get_v_master_cells
Definition: get_master.f:170
bc::bc_area
double precision, dimension(dimension_bc) bc_area
Definition: bc_mod.f:245
reasssign_quadric
subroutine reasssign_quadric(x1, x2, x3, GROUP, Q_ID)
Definition: define_quadrics.f:621
cutcell::unit_cut_cell_log
integer unit_cut_cell_log
Definition: cutcell_mod.f:8
cutcell::snap
logical, dimension(:), allocatable snap
Definition: cutcell_mod.f:413