File: N:\mfix\model\cartesian_grid\get_cut_cell_volume_area.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE GET_CUT_CELL_VOLUME_AND_AREAS(IJK,TYPE_OF_CELL,N_NODES,CONNECT,X_NP,Y_NP,Z_NP)
16
17 USE bc
18 USE compar, ONLY: ijkend3, istart1, jstart1, kstart1, mype
19 USE cut_cell_preproc, only: eval_f
20 USE cutcell
21 USE exit, only: mfix_exit
22 USE functions, only: funijk
23 USE geometry, ONLY: do_k, no_k, ayz, axz, axy, ayz_u, axy_u, axz_u, ayz_v, axz_v, axy_v, ayz_w, axz_w, axy_w, vol, vol_u, vol_v, vol_w, dx, dy, dz, zlength, flag, flag_n, flag_t
24 USE indices, only: i_of, j_of, k_of
25 USE param, only: dimension_3
26 USE param1, only: undefined
27 USE param1, only: zero
28 USE polygon, ONLY: n_polygon
29 USE quadric, ONLY: n_quadric, tol_f, bc_id_q, dim_quadric
30 USE stl, ONLY: n_facet_at, list_facet_at, bc_id_stl_face
31
32 IMPLICIT NONE
33 CHARACTER (LEN=*) :: TYPE_OF_CELL
34 INTEGER :: I,J,K,IJK,L,NODE
35 INTEGER :: IM,JM,KM,IMJK,IJMK,IJKM
36 INTEGER :: NN,N_N1,N_N2,Q_ID
37 INTEGER :: N_CUT_FACE_NODES
38 INTEGER :: N_EAST_FACE_NODES,N_NORTH_FACE_NODES,N_TOP_FACE_NODES
39 INTEGER :: N_WEST_FACE_NODES,N_SOUTH_FACE_NODES,N_BOTTOM_FACE_NODES
40 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_CUT_FACE_NODES
41 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_EAST_FACE_NODES
42 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_NORTH_FACE_NODES
43 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_TOP_FACE_NODES
44 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_WEST_FACE_NODES
45 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_SOUTH_FACE_NODES
46 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_BOTTOM_FACE_NODES
47 DOUBLE PRECISION :: X_COPY,Y_COPY,Z_COPY,F_COPY,F_Q,F_MIN
48 DOUBLE PRECISION :: X_MEAN,Y_MEAN,Z_MEAN
49 DOUBLE PRECISION :: AREA_EAST,AREA_NORTH,AREA_TOP
50 DOUBLE PRECISION :: AREA_WEST,AREA_SOUTH,AREA_BOTTOM
51 DOUBLE PRECISION :: CUT_AREA
52 DOUBLE PRECISION, DIMENSION(3) :: CENTROID_EAST,CENTROID_NORTH,CENTROID_TOP
53 DOUBLE PRECISION, DIMENSION(3) :: CENTROID_WEST,CENTROID_SOUTH,CENTROID_BOTTOM
54 DOUBLE PRECISION, DIMENSION(3) :: CENTROID_CUT
55 DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz,VOLUME,NORM
56 DOUBLE PRECISION :: DIFFERENCE
57 INTEGER :: N_NODES
58 INTEGER :: NDIFF,SPACE
59 INTEGER, DIMENSION(DIMENSION_3,15) :: CONNECT
60 DOUBLE PRECISION, DIMENSION(DIMENSION_MAX_CUT_CELL) :: X_NP,Y_NP,Z_NP
61
62 DOUBLE PRECISION, DIMENSION(DIM_QUADRIC) :: F_CUT_FACE
63
64 INTEGER, DIMENSION(15) :: TEMP_CONNECTIVITY
65
66 LOGICAL :: CLIP_FLAG,INSIDE_FACET
67 INTEGER :: BCID,N_BC,QID_FMIN,BCID2,NF
68
69 LOGICAL :: CORNER_POINT
70 INTEGER :: NODE_OF_CORNER,IERROR
71
72
73
74
75
76
77 = I_OF(IJK)
78 J = J_OF(IJK)
79 K = K_OF(IJK)
80
81 IM = I - 1
82 JM = J - 1
83 KM = K - 1
84
85 IMJK = FUNIJK(IM,J,K)
86 IJMK = FUNIJK(I,JM,K)
87 IJKM = FUNIJK(I,J,KM)
88
89 N_EAST_FACE_NODES = 0
90 N_NORTH_FACE_NODES = 0
91 N_TOP_FACE_NODES = 0
92
93 N_CUT_FACE_NODES = 0
94
95 N_WEST_FACE_NODES = 0
96 N_SOUTH_FACE_NODES = 0
97 N_BOTTOM_FACE_NODES = 0
98
99 X_MEAN = ZERO
100 Y_MEAN = ZERO
101 Z_MEAN = ZERO
102
103 IF(NO_K) THEN
104 N_N1 = 5
105 N_N2 = 8
106 ELSE
107 N_N1 = 1
108 N_N2 = 8
109 ENDIF
110
111 DO L = 1, N_NODES
112 IF(CONNECT(IJK,L)>IJKEND3) THEN
113 = X_NP(CONNECT(IJK,L)-IJKEND3)
114 Y_COPY = Y_NP(CONNECT(IJK,L)-IJKEND3)
115 Z_COPY = Z_NP(CONNECT(IJK,L)-IJKEND3)
116 CORNER_POINT = .FALSE.
117
118 ELSE
119 = .TRUE.
120 DO NODE = N_N1,N_N2
121 IF(CONNECT(IJK,L) == IJK_OF_NODE(NODE)) THEN
122 NODE_OF_CORNER = NODE
123 X_COPY = X_NODE(NODE)
124 Y_COPY = Y_NODE(NODE)
125 Z_COPY = Z_NODE(NODE)
126
127 IF(SNAP(IJK_OF_NODE(NODE))) THEN
128 N_CUT_FACE_NODES = N_CUT_FACE_NODES + 1
129 COORD_CUT_FACE_NODES(1,N_CUT_FACE_NODES) = X_COPY
130 COORD_CUT_FACE_NODES(2,N_CUT_FACE_NODES) = Y_COPY
131 COORD_CUT_FACE_NODES(3,N_CUT_FACE_NODES) = Z_COPY
132
133 IF(N_QUADRIC>0) THEN
134 F_MIN = UNDEFINED
135 QID_FMIN = 0
136 DO Q_ID = 1, N_QUADRIC
137 CALL GET_F_QUADRIC(X_COPY,Y_COPY,Z_COPY,Q_ID,F_Q,CLIP_FLAG)
138 IF(DABS(F_Q) < F_MIN) THEN
139 F_MIN = DABS(F_Q)
140 QID_FMIN = Q_ID
141 ENDIF
142 ENDDO
143 BC_ID(IJK) = BC_ID_Q(QID_FMIN)
144 ENDIF
145
146 ENDIF
147
148 ENDIF
149 END DO
150 ENDIF
151
152 X_MEAN = X_MEAN + X_COPY
153 Y_MEAN = Y_MEAN + Y_COPY
154 Z_MEAN = Z_MEAN + Z_COPY
155
156 IF(CORNER_POINT) THEN
157 Q_ID = 1
158 CALL EVAL_F('QUADRIC',X_COPY,Y_COPY,Z_COPY,Q_ID,F_COPY,CLIP_FLAG)
159
160 CALL EVAL_F('POLYGON',X_COPY,Y_COPY,Z_COPY,N_POLYGON,F_COPY,CLIP_FLAG)
161
162 CALL EVAL_F('USR_DEF',X_COPY,Y_COPY,Z_COPY,N_USR_DEF,F_COPY,CLIP_FLAG)
163 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE_OF_CORNER,F_COPY,CLIP_FLAG,BCID2)
164 ELSE
165 F_COPY = ZERO
166 ENDIF
167
168 IF (DABS(F_COPY) < TOL_F ) THEN
169 = N_CUT_FACE_NODES + 1
170 COORD_CUT_FACE_NODES(1,N_CUT_FACE_NODES) = X_COPY
171 COORD_CUT_FACE_NODES(2,N_CUT_FACE_NODES) = Y_COPY
172 COORD_CUT_FACE_NODES(3,N_CUT_FACE_NODES) = Z_COPY
173 ENDIF
174
175 IF (X_COPY > (X_NODE(8) - TOL_F) ) THEN
176 = N_EAST_FACE_NODES + 1
177 COORD_EAST_FACE_NODES(1,N_EAST_FACE_NODES) = X_COPY
178 COORD_EAST_FACE_NODES(2,N_EAST_FACE_NODES) = Y_COPY
179 COORD_EAST_FACE_NODES(3,N_EAST_FACE_NODES) = Z_COPY
180 ENDIF
181
182 IF (Y_COPY > (Y_NODE(8) - TOL_F) ) THEN
183 = N_NORTH_FACE_NODES + 1
184 COORD_NORTH_FACE_NODES(1,N_NORTH_FACE_NODES) = X_COPY
185 COORD_NORTH_FACE_NODES(2,N_NORTH_FACE_NODES) = Y_COPY
186 COORD_NORTH_FACE_NODES(3,N_NORTH_FACE_NODES) = Z_COPY
187 ENDIF
188
189 IF (Z_COPY > (Z_NODE(8) - TOL_F) ) THEN
190 = N_TOP_FACE_NODES + 1
191 COORD_TOP_FACE_NODES(1,N_TOP_FACE_NODES) = X_COPY
192 COORD_TOP_FACE_NODES(2,N_TOP_FACE_NODES) = Y_COPY
193 COORD_TOP_FACE_NODES(3,N_TOP_FACE_NODES) = Z_COPY
194 ENDIF
195
196 IF(I>=ISTART1) THEN
197 IF (X_COPY < (X_NODE(1) + TOL_F) ) THEN
198 = N_WEST_FACE_NODES + 1
199 COORD_WEST_FACE_NODES(1,N_WEST_FACE_NODES) = X_COPY
200 COORD_WEST_FACE_NODES(2,N_WEST_FACE_NODES) = Y_COPY
201 COORD_WEST_FACE_NODES(3,N_WEST_FACE_NODES) = Z_COPY
202 ENDIF
203 ENDIF
204
205 IF(J>=JSTART1) THEN
206 IF (Y_COPY < (Y_NODE(1) + TOL_F) ) THEN
207 = N_SOUTH_FACE_NODES + 1
208 COORD_SOUTH_FACE_NODES(1,N_SOUTH_FACE_NODES) = X_COPY
209 COORD_SOUTH_FACE_NODES(2,N_SOUTH_FACE_NODES) = Y_COPY
210 COORD_SOUTH_FACE_NODES(3,N_SOUTH_FACE_NODES) = Z_COPY
211 ENDIF
212 ENDIF
213
214 IF(DO_K) THEN
215 IF(K>=KSTART1) THEN
216 IF (Z_COPY < (Z_NODE(1) + TOL_F) ) THEN
217 = N_BOTTOM_FACE_NODES + 1
218 COORD_BOTTOM_FACE_NODES(1,N_BOTTOM_FACE_NODES) = X_COPY
219 COORD_BOTTOM_FACE_NODES(2,N_BOTTOM_FACE_NODES) = Y_COPY
220 COORD_BOTTOM_FACE_NODES(3,N_BOTTOM_FACE_NODES) = Z_COPY
221 ENDIF
222 ENDIF
223 ENDIF
224
225 END DO
226
227 X_MEAN = X_MEAN / N_NODES
228 Y_MEAN = Y_MEAN / N_NODES
229 Z_MEAN = Z_MEAN / N_NODES
230
231
232 IF( N_CUT_FACE_NODES == N_EAST_FACE_NODES) THEN
233 DIFFERENCE = ZERO
234 DO NDIFF = 1,N_CUT_FACE_NODES
235 DO SPACE = 1,3
236 DIFFERENCE = DIFFERENCE +DABS(COORD_CUT_FACE_NODES(SPACE,NDIFF) &
237 - COORD_EAST_FACE_NODES(SPACE,NDIFF))
238 END DO
239 END DO
240 IF( DIFFERENCE < TOL_F ) THEN
241 IF(PRINT_WARNINGS) THEN
242 WRITE(*,*)'WARNING: CUT FACE IS IDENTICAL TO EAST FACE AT IJK = ',IJK
243 WRITE(*,*)'TYPE OF CELL = ',TYPE_OF_CELL
244 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
245 WRITE(*,*)'REMOVING EAST FACE.'
246 ENDIF
247 N_EAST_FACE_NODES = 0
248 ENDIF
249 ENDIF
250
251 IF( N_CUT_FACE_NODES == N_WEST_FACE_NODES) THEN
252 DIFFERENCE = ZERO
253 DO NDIFF = 1,N_CUT_FACE_NODES
254 DO SPACE = 1,3
255 DIFFERENCE = DIFFERENCE +DABS(COORD_CUT_FACE_NODES(SPACE,NDIFF) &
256 - COORD_WEST_FACE_NODES(SPACE,NDIFF))
257 END DO
258 END DO
259 IF( DIFFERENCE < TOL_F ) THEN
260 IF(PRINT_WARNINGS) THEN
261 WRITE(*,*)'WARNING: CUT FACE IS IDENTICAL TO WEST FACE AT IJK = ',IJK
262 WRITE(*,*)'TYPE OF CELL = ',TYPE_OF_CELL
263 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
264 WRITE(*,*)'REMOVING WEST FACE.'
265 ENDIF
266 N_WEST_FACE_NODES = 0
267 ENDIF
268 ENDIF
269
270 IF( N_CUT_FACE_NODES == N_NORTH_FACE_NODES) THEN
271 DIFFERENCE = ZERO
272 DO NDIFF = 1,N_CUT_FACE_NODES
273 DO SPACE = 1,3
274 DIFFERENCE = DIFFERENCE +DABS(COORD_CUT_FACE_NODES(SPACE,NDIFF) &
275 - COORD_NORTH_FACE_NODES(SPACE,NDIFF))
276 END DO
277 END DO
278 IF( DIFFERENCE < TOL_F ) THEN
279 IF(PRINT_WARNINGS) THEN
280 WRITE(*,*)'WARNING: CUT FACE IS IDENTICAL TO NORTH FACE AT IJK = ',IJK
281 WRITE(*,*)'TYPE OF CELL = ',TYPE_OF_CELL
282 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
283 WRITE(*,*)'REMOVING NORTH FACE.'
284 ENDIF
285 N_NORTH_FACE_NODES = 0
286 ENDIF
287 ENDIF
288
289 IF( N_CUT_FACE_NODES == N_SOUTH_FACE_NODES) THEN
290 DIFFERENCE = ZERO
291 DO NDIFF = 1,N_CUT_FACE_NODES
292 DO SPACE = 1,3
293 DIFFERENCE = DIFFERENCE +DABS(COORD_CUT_FACE_NODES(SPACE,NDIFF) &
294 - COORD_SOUTH_FACE_NODES(SPACE,NDIFF))
295 END DO
296 END DO
297 IF( DIFFERENCE < TOL_F ) THEN
298 IF(PRINT_WARNINGS) THEN
299 WRITE(*,*)'WARNING: CUT FACE IS IDENTICAL TO SOUTH FACE AT IJK = ',IJK
300 WRITE(*,*)'TYPE OF CELL = ',TYPE_OF_CELL
301 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
302 WRITE(*,*)'REMOVING SOUTH FACE.'
303 ENDIF
304 N_SOUTH_FACE_NODES = 0
305 ENDIF
306 ENDIF
307
308 IF((NO_K).AND.(N_TOP_FACE_NODES <=2)) THEN
309 IF(PRINT_WARNINGS) THEN
310 WRITE(*,*)'WARNING: TOP FACE HAS ONLY TWO NODES AT IJK = ',IJK
311 WRITE(*,*)'TYPE OF CELL = ',TYPE_OF_CELL
312 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
313 WRITE(*,*)'REMOVING CUT CELL.'
314 ENDIF
315
316 SELECT CASE (TYPE_OF_CELL)
317 CASE('SCALAR')
318 BLOCKED_CELL_AT(IJK) = .TRUE.
319 CASE('U_MOMENTUM')
320 BLOCKED_U_CELL_AT(IJK) = .TRUE.
321 CASE('V_MOMENTUM')
322 BLOCKED_V_CELL_AT(IJK) = .TRUE.
323 CASE('W_MOMENTUM')
324 BLOCKED_W_CELL_AT(IJK) = .TRUE.
325 CASE DEFAULT
326 WRITE(*,*)'SUBROUTINE: GET_CUT_CELL_VOLUME_AND_AREAS'
327 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
328 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
329 WRITE(*,*)'SCALAR'
330 WRITE(*,*)'U_MOMENTUM'
331 WRITE(*,*)'V_MOMENTUM'
332 WRITE(*,*)'W_MOMENTUM'
333 CALL MFIX_EXIT(myPE)
334 END SELECT
335 RETURN
336
337 ELSE
338 IF( N_CUT_FACE_NODES == N_TOP_FACE_NODES) THEN
339 DIFFERENCE = ZERO
340 DO NDIFF = 1,N_CUT_FACE_NODES
341 DO SPACE = 1,3
342 DIFFERENCE = DIFFERENCE +DABS(COORD_CUT_FACE_NODES(SPACE,NDIFF) &
343 - COORD_TOP_FACE_NODES(SPACE,NDIFF))
344 END DO
345 END DO
346 IF( DIFFERENCE < TOL_F ) THEN
347 IF(PRINT_WARNINGS) THEN
348 WRITE(*,*)'WARNING: CUT FACE IS IDENTICAL TO TOP FACE AT IJK = ',IJK
349 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
350 WRITE(*,*)'REMOVING TOP FACE.'
351 ENDIF
352 N_TOP_FACE_NODES = 0
353 ENDIF
354 ENDIF
355
356 IF( N_CUT_FACE_NODES == N_BOTTOM_FACE_NODES) THEN
357 DIFFERENCE = ZERO
358 DO NDIFF = 1,N_CUT_FACE_NODES
359 DO SPACE = 1,3
360 DIFFERENCE = DIFFERENCE +DABS(COORD_CUT_FACE_NODES(SPACE,NDIFF) &
361 - COORD_BOTTOM_FACE_NODES(SPACE,NDIFF))
362 END DO
363 END DO
364 IF( DIFFERENCE < TOL_F ) THEN
365 IF(PRINT_WARNINGS) THEN
366 WRITE(*,*)'WARNING: CUT FACE IS IDENTICAL TO BOTTOM FACE AT IJK = ',IJK
367 WRITE(*,*)'THIS USUALLY OCCURS WHEN QUADRIC SURFACE IS ALIGNED WITH GRID.'
368 WRITE(*,*)'REMOVING BOTTOM FACE.'
369 ENDIF
370 N_BOTTOM_FACE_NODES = 0
371 ENDIF
372 ENDIF
373 ENDIF
374
375 IERROR = 0
376
377 CALL GET_POLYGON_AREA_AND_CENTROID(N_EAST_FACE_NODES,COORD_EAST_FACE_NODES,AREA_EAST,CENTROID_EAST,IERROR)
378
379 CALL GET_POLYGON_AREA_AND_CENTROID(N_NORTH_FACE_NODES,COORD_NORTH_FACE_NODES,AREA_NORTH,CENTROID_NORTH,IERROR)
380
381 CALL GET_POLYGON_AREA_AND_CENTROID(N_TOP_FACE_NODES,COORD_TOP_FACE_NODES,AREA_TOP,CENTROID_TOP,IERROR)
382
383 CALL GET_POLYGON_AREA_AND_CENTROID(N_CUT_FACE_NODES,COORD_CUT_FACE_NODES,CUT_AREA,CENTROID_CUT,IERROR)
384
385 CALL STORE_CUT_FACE_INFO(IJK,TYPE_OF_CELL,N_CUT_FACE_NODES,COORD_CUT_FACE_NODES,X_MEAN,Y_MEAN,Z_MEAN)
386
387 CALL GET_DEL_H(IJK,TYPE_OF_CELL,X_MEAN,Y_MEAN,Z_MEAN,Del_H,Nx,Ny,Nz)
388
389 IF(DEL_H == UNDEFINED) DEL_H = ZERO
390
391 IF(I==ISTART1) CALL GET_POLYGON_AREA_AND_CENTROID(N_WEST_FACE_NODES,COORD_WEST_FACE_NODES,AREA_WEST,CENTROID_WEST,IERROR)
392
393 IF(J==JSTART1) CALL GET_POLYGON_AREA_AND_CENTROID(N_SOUTH_FACE_NODES,COORD_SOUTH_FACE_NODES,AREA_SOUTH,CENTROID_SOUTH,IERROR)
394
395 IF(DO_K.AND.K==KSTART1) CALL GET_POLYGON_AREA_AND_CENTROID(N_BOTTOM_FACE_NODES,COORD_BOTTOM_FACE_NODES,AREA_BOTTOM,&
396 CENTROID_BOTTOM,IERROR)
397
398 SELECT CASE (TYPE_OF_CELL)
399 CASE('SCALAR')
400 IF(I>ISTART1) AREA_WEST = AYZ(IMJK)
401 IF(J>JSTART1) AREA_SOUTH = AXZ(IJMK)
402 IF(K>KSTART1) AREA_BOTTOM = AXY(IJKM)
403
404 AYZ(IJK) = AREA_EAST
405 AXZ(IJK) = AREA_NORTH
406 AXY(IJK) = AREA_TOP
407 Area_CUT(IJK) = CUT_AREA
408
409
410 (IJK,1) = AREA_EAST - AREA_WEST
411 NORMAL_S(IJK,2) = AREA_NORTH - AREA_SOUTH
412
413 IF(DO_K) THEN
414 NORMAL_S(IJK,3) = AREA_TOP - AREA_BOTTOM
415 ELSE
416 NORMAL_S(IJK,3) = ZERO
417 ENDIF
418
419 NORM = sqrt(NORMAL_S(IJK,1)**2 + NORMAL_S(IJK,2)**2 + NORMAL_S(IJK,3)**2)
420
421 IF(NORM==ZERO) NORM = UNDEFINED
422
423 NORMAL_S(IJK,:) = NORMAL_S(IJK,:) / NORM
424
425 REFP_S(IJK,:) = CENTROID_CUT(:)
426
427 CALL GET_DEL_H(IJK,TYPE_OF_CELL,X_MEAN,Y_MEAN,Z_MEAN,Del_H,Nx,Ny,Nz)
428 IF(DEL_H == UNDEFINED) DEL_H = ZERO
429 DELH_Scalar(IJK) = DEL_H
430
431 IF(NO_K) THEN
432 VOL(IJK) = AXY(IJK) * ZLENGTH
433 ELSE
434
435
436
437 = AREA_EAST * (X_NODE(8) - X_MEAN) &
438 + AREA_NORTH * (Y_NODE(8) - Y_MEAN) &
439 + AREA_TOP * (Z_NODE(8) - Z_MEAN) &
440 + AREA_WEST * (X_MEAN - X_NODE(1)) &
441 + AREA_SOUTH * (Y_MEAN - Y_NODE(1)) &
442 + AREA_BOTTOM * (Z_MEAN - Z_NODE(1)) &
443 + CUT_AREA * DEL_H
444
445 VOL(IJK) = VOLUME / 3.0D0
446
447 ENDIF
448
449 IF(IERROR==1) print*,'scalar IERROR',IJK
450 IF(IERROR==1) VOL(IJK) = ZERO
451
452 IF(VOL(IJK)==ZERO) THEN
453
454 IF(PRINT_WARNINGS) THEN
455 PRINT*,'WARNING: ZERO VOLUME CUT CELL DETECTED AT IJK =',IJK
456 PRINT*,'REMOVING CUT CELL FROM COMPUTATION (FLAGGED AS BLOCKED CELL)'
457 ENDIF
458
459 FLAG(IJK) = 100
460 CUT_CELL_AT(IJK) = .FALSE.
461 BLOCKED_CELL_AT(IJK) = .TRUE.
462 STANDARD_CELL_AT(IJK) = .FALSE.
463 AYZ(IJK) = ZERO
464 AXZ(IJK) = ZERO
465 AXY(IJK) = ZERO
466 VOL(IJK) = ZERO
467
468 ELSEIF(VOL(IJK) < TOL_SMALL_CELL * DX(I)*DY(J)*DZ(K) ) THEN
469 SMALL_CELL_AT(IJK) = .TRUE.
470 NUMBER_OF_SMALL_CELLS = NUMBER_OF_SMALL_CELLS + 1
471
472
473 IF(PRINT_WARNINGS) THEN
474 PRINT*,'WARNING: SMALL SCALAR CELL CUT CELL DETECTED AT IJK =',IJK
475 PRINT*,'REMOVING CUT CELL FROM COMPUTATION (FLAGGED AS BLOCKED CELL)'
476 ENDIF
477
478 BLOCKED_CELL_AT(IJK) = .TRUE.
479
480 ENDIF
481
482 IF(AREA_EAST > TOL_SMALL_AREA*DY(J)*DZ(K)) THEN
483 WALL_U_AT(IJK) = .FALSE.
484 X_U(IJK) = CENTROID_EAST(1)
485 Y_U(IJK) = CENTROID_EAST(2)
486 Z_U(IJK) = CENTROID_EAST(3)
487 ENDIF
488
489 IF(AREA_NORTH > TOL_SMALL_AREA*DX(I)*DZ(K)) THEN
490 WALL_V_AT(IJK) = .FALSE.
491 X_V(IJK) = CENTROID_NORTH(1)
492 Y_V(IJK) = CENTROID_NORTH(2)
493 Z_V(IJK) = CENTROID_NORTH(3)
494 ENDIF
495
496 IF(AREA_TOP > TOL_SMALL_AREA*DX(I)*DY(J)) THEN
497 WALL_W_AT(IJK) = .FALSE.
498 X_W(IJK) = CENTROID_TOP(1)
499 Y_W(IJK) = CENTROID_TOP(2)
500 Z_W(IJK) = CENTROID_TOP(3)
501 ENDIF
502
503 IF(N_QUADRIC>0) THEN
504 F_CUT_FACE = UNDEFINED
505 N_BC = 0
506
507 DO Q_ID = 1, N_QUADRIC
508 DO NODE = 1,N_CUT_FACE_NODES
509
510 CALL GET_F_QUADRIC(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
511 Q_ID,F_CUT_FACE(Q_ID),CLIP_FLAG)
512
513 IF(DABS(F_CUT_FACE(Q_ID)) < TOL_F) THEN
514 IF(BC_ID(IJK)/=BC_ID_Q(Q_ID)) N_BC = N_BC + 1
515 BC_ID(IJK) = BC_ID_Q(Q_ID)
516 ENDIF
517 ENDDO
518 IF(BC_ID(IJK)>0) THEN
519
520 IF(BC_TYPE_ENUM(BC_ID(IJK)) == CG_NSW) EXIT
521 ENDIF
522 ENDDO
523 ENDIF
524
525 IF(N_POLYGON>0) THEN
526 DO NODE = 1,N_CUT_FACE_NODES
527 CALL EVAL_POLY_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
528 Q_ID,F_CUT_FACE(1),CLIP_FLAG,BCID)
529 IF(F_CUT_FACE(1)==ZERO) THEN
530 BC_ID(IJK) = BCID
531 IF(BC_ID(IJK)>0) THEN
532 IF(BC_TYPE_ENUM(BC_ID(IJK)) == CG_MI) EXIT
533 ENDIF
534 ENDIF
535 ENDDO
536 ENDIF
537
538 IF(N_USR_DEF>0) THEN
539 DO NODE = 1,N_CUT_FACE_NODES
540 CALL EVAL_USR_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
541 BCID,F_CUT_FACE(1),CLIP_FLAG)
542 IF(DABS(F_CUT_FACE(1)) < TOL_F) THEN
543 BC_ID(IJK) = BCID
544 IF(BC_ID(IJK)>0) THEN
545 IF(BC_TYPE_ENUM(BC_ID(IJK)) == CG_MI) EXIT
546 ENDIF
547 ENDIF
548 ENDDO
549 ENDIF
550
551 IF(USE_STL.OR.USE_MSH) THEN
552 DO NODE = 1,N_CUT_FACE_NODES
553 DO NN=1,N_FACET_AT(IJK)
554 NF=LIST_FACET_AT(IJK,NN)
555 CALL IS_POINT_INSIDE_FACET(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),&
556 COORD_CUT_FACE_NODES(3,NODE),NF,INSIDE_FACET)
557 IF(INSIDE_FACET) THEN
558 BC_ID(IJK) = BC_ID_STL_FACE(NF)
559 IF(BC_ID(IJK)>0) THEN
560 IF(BC_TYPE_ENUM(BC_ID(IJK)) == CG_MI) EXIT
561 ENDIF
562 ENDIF
563 ENDDO
564 ENDDO
565 ENDIF
566
567 IF (BC_ID(IJK)>0) THEN
568
569 IF (.not. IS_CG(BC_TYPE_ENUM(BC_ID(IJK)))) THEN
570 WRITE(*,'(/1X,A)')'ERROR: INVALID BOUNDARY CONDITION ASSIGNED TO A CUT CELL.'
571 WRITE(*,'(/1X,A,I9,I6)')'IJK, MyPE = ', IJK,myPE
572 WRITE(*,'(1X,A,I6)')'BC_ID = ',BC_ID(IJK)
573 WRITE(*,'(1X,A,A)')'BC_TYPE = ',BC_TYPE_ENUM(BC_ID(IJK))
574 WRITE(*,'(1X,A)')'VALID BC_TYPE FOR CUT CELLS ARE: CG_NSW, CG_FSW, CG_PSW, CG_MI, and CG_PO'
575 WRITE(*,'(/1X,A)')'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
576 call mfix_exit(myPE)
577 ENDIF
578
579 ENDIF
580
581
582
583 IF(NO_K) THEN
584 DO NN = 1,N_TOP_FACE_NODES
585 TEMP_CONNECTIVITY(NN) = CONNECT(IJK,ORDER(NN))
586 END DO
587 CONNECT(IJK,:) = TEMP_CONNECTIVITY
588 ENDIF
589
590 CASE('U_MOMENTUM')
591
592 IF(I>ISTART1) AREA_WEST = AYZ_U(IMJK)
593 IF(J>JSTART1) AREA_SOUTH = AXZ_U(IJMK)
594 IF(K>KSTART1) AREA_BOTTOM = AXY_U(IJKM)
595
596 AYZ_U(IJK) = AREA_EAST
597 AXZ_U(IJK) = AREA_NORTH
598 AXY_U(IJK) = AREA_TOP
599 Area_U_CUT(IJK) = CUT_AREA
600
601 IF(AREA_EAST > ZERO) THEN
602 X_U_ec(IJK) = CENTROID_EAST(1)
603 Y_U_ec(IJK) = CENTROID_EAST(2)
604 Z_U_ec(IJK) = CENTROID_EAST(3)
605 ENDIF
606
607 IF(AREA_NORTH > ZERO) THEN
608 X_U_nc(IJK) = CENTROID_NORTH(1)
609 Y_U_nc(IJK) = CENTROID_NORTH(2)
610 Z_U_nc(IJK) = CENTROID_NORTH(3)
611 ENDIF
612
613 IF(AREA_TOP > ZERO) THEN
614 X_U_tc(IJK) = CENTROID_TOP(1)
615 Y_U_tc(IJK) = CENTROID_TOP(2)
616 Z_U_tc(IJK) = CENTROID_TOP(3)
617 ENDIF
618
619
620
621 (IJK,1) = AREA_EAST - AREA_WEST
622 NORMAL_U(IJK,2) = AREA_NORTH - AREA_SOUTH
623
624 IF(DO_K) THEN
625 NORMAL_U(IJK,3) = AREA_TOP - AREA_BOTTOM
626 ELSE
627 NORMAL_U(IJK,3) = ZERO
628 ENDIF
629
630 NORM = sqrt(NORMAL_U(IJK,1)**2 + NORMAL_U(IJK,2)**2 + NORMAL_U(IJK,3)**2)
631
632 IF(NORM==ZERO) NORM = UNDEFINED
633
634 NORMAL_U(IJK,:) = NORMAL_U(IJK,:) / NORM
635
636 REFP_U(IJK,:) = CENTROID_CUT(:)
637
638 CALL GET_DEL_H(IJK,TYPE_OF_CELL,X_MEAN,Y_MEAN,Z_MEAN,Del_H,Nx,Ny,Nz)
639 IF(DEL_H == UNDEFINED) DEL_H = ZERO
640
641 IF(NO_K) THEN
642 VOL_U(IJK) = AXY_U(IJK) * ZLENGTH
643 ELSE
644
645
646
647 = AREA_EAST * (X_NODE(8) - X_MEAN) &
648 + AREA_NORTH * (Y_NODE(8) - Y_MEAN) &
649 + AREA_TOP * (Z_NODE(8) - Z_MEAN) &
650 + AREA_WEST * (X_MEAN - X_NODE(1)) &
651 + AREA_SOUTH * (Y_MEAN - Y_NODE(1)) &
652 + AREA_BOTTOM * (Z_MEAN - Z_NODE(1)) &
653 + CUT_AREA * DEL_H
654
655 VOL_U(IJK) = VOLUME / 3.0D0
656 ENDIF
657
658 IF(IERROR==1) print*,'U_VEL IERROR',IJK
659 IF(IERROR==1) AYZ(IJK) = ZERO
660
661 IF(AYZ(IJK) <= TOL_SMALL_AREA*DY(J)*DZ(K)) THEN
662 WALL_U_AT(IJK) = .TRUE.
663
664 NUMBER_OF_U_WALL_CELLS = NUMBER_OF_U_WALL_CELLS + 1
665
666 X_U(IJK) = CENTROID_CUT(1)
667 Y_U(IJK) = CENTROID_CUT(2)
668 Z_U(IJK) = CENTROID_CUT(3)
669 ENDIF
670
671
672 IF(N_QUADRIC>0) THEN
673 F_CUT_FACE = UNDEFINED
674 N_BC = 0
675 BC_U_ID(IJK) = 0
676 DO Q_ID = 1, N_QUADRIC
677 DO NODE = 1,N_CUT_FACE_NODES
678
679 CALL GET_F_QUADRIC(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
680 Q_ID,F_CUT_FACE(Q_ID),CLIP_FLAG)
681
682 IF(DABS(F_CUT_FACE(Q_ID)) < TOL_F) THEN
683 IF(BC_U_ID(IJK)/=BC_ID_Q(Q_ID)) N_BC = N_BC + 1
684 BC_U_ID(IJK) = BC_ID_Q(Q_ID)
685 ENDIF
686 ENDDO
687 IF(BC_U_ID(IJK)>0) THEN
688 IF(BC_TYPE_ENUM(BC_U_ID(IJK)) == CG_MI) EXIT
689
690 ENDIF
691
692 ENDDO
693
694 ENDIF
695
696 IF(N_POLYGON>0) THEN
697 DO NODE = 1,N_CUT_FACE_NODES
698 CALL EVAL_POLY_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
699 Q_ID,F_CUT_FACE(1),CLIP_FLAG,BCID)
700 IF(F_CUT_FACE(1)==ZERO) THEN
701 BC_U_ID(IJK) = BCID
702 IF(BC_U_ID(IJK)>0) THEN
703 IF(BC_TYPE_ENUM(BC_U_ID(IJK)) == CG_MI) EXIT
704 ENDIF
705 ENDIF
706 ENDDO
707 ENDIF
708
709
710 IF(N_USR_DEF>0) THEN
711 DO NODE = 1,N_CUT_FACE_NODES
712 CALL EVAL_USR_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
713 BCID,F_CUT_FACE(1),CLIP_FLAG)
714 IF(DABS(F_CUT_FACE(1)) < TOL_F) THEN
715 BC_U_ID(IJK) = BCID
716 IF(BC_U_ID(IJK)>0) THEN
717 IF(BC_TYPE_ENUM(BC_U_ID(IJK)) == CG_MI) EXIT
718 ENDIF
719 ENDIF
720 ENDDO
721 ENDIF
722
723
724 IF(USE_STL.OR.USE_MSH) THEN
725 DO NODE = 1,N_CUT_FACE_NODES
726 DO NN=1,N_FACET_AT(IJK)
727 NF=LIST_FACET_AT(IJK,NN)
728 CALL IS_POINT_INSIDE_FACET(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),&
729 COORD_CUT_FACE_NODES(3,NODE),NF,INSIDE_FACET)
730 IF(INSIDE_FACET) THEN
731 BC_U_ID(IJK) = BC_ID_STL_FACE(NF)
732 IF(BC_U_ID(IJK)>0) THEN
733 IF(BC_TYPE_ENUM(BC_U_ID(IJK)) == CG_MI) EXIT
734 ENDIF
735 ENDIF
736 ENDDO
737 ENDDO
738 ENDIF
739
740 CASE('V_MOMENTUM')
741
742 IF(I>ISTART1) AREA_WEST = AYZ_V(IMJK)
743 IF(J>JSTART1) AREA_SOUTH = AXZ_V(IJMK)
744 IF(K>KSTART1) AREA_BOTTOM = AXY_V(IJKM)
745
746 AYZ_V(IJK) = AREA_EAST
747 AXZ_V(IJK) = AREA_NORTH
748 AXY_V(IJK) = AREA_TOP
749 Area_V_CUT(IJK) = CUT_AREA
750
751 IF(AREA_EAST > ZERO) THEN
752 X_V_ec(IJK) = CENTROID_EAST(1)
753 Y_V_ec(IJK) = CENTROID_EAST(2)
754 Z_V_ec(IJK) = CENTROID_EAST(3)
755 ENDIF
756
757 IF(AREA_NORTH > ZERO) THEN
758 X_V_nc(IJK) = CENTROID_NORTH(1)
759 Y_V_nc(IJK) = CENTROID_NORTH(2)
760 Z_V_nc(IJK) = CENTROID_NORTH(3)
761 ENDIF
762
763 IF(AREA_TOP > ZERO) THEN
764 X_V_tc(IJK) = CENTROID_TOP(1)
765 Y_V_tc(IJK) = CENTROID_TOP(2)
766 Z_V_tc(IJK) = CENTROID_TOP(3)
767 ENDIF
768
769
770 (IJK,1) = AREA_EAST - AREA_WEST
771 NORMAL_V(IJK,2) = AREA_NORTH - AREA_SOUTH
772
773 IF(DO_K) THEN
774 NORMAL_V(IJK,3) = AREA_TOP - AREA_BOTTOM
775 ELSE
776 NORMAL_V(IJK,3) = ZERO
777 ENDIF
778
779
780 NORM = sqrt(NORMAL_V(IJK,1)**2 + NORMAL_V(IJK,2)**2 + NORMAL_V(IJK,3)**2)
781
782 IF(NORM==ZERO) NORM = UNDEFINED
783
784 NORMAL_V(IJK,:) = NORMAL_V(IJK,:) / NORM
785
786 REFP_V(IJK,:) = CENTROID_CUT(:)
787
788 CALL GET_DEL_H(IJK,TYPE_OF_CELL,X_MEAN,Y_MEAN,Z_MEAN,Del_H,Nx,Ny,Nz)
789 IF(DEL_H == UNDEFINED) DEL_H = ZERO
790
791 IF(NO_K) THEN
792 VOL_V(IJK) = AXY_V(IJK) * ZLENGTH
793 ELSE
794
795
796
797 = AREA_EAST * (X_NODE(8) - X_MEAN) &
798 + AREA_NORTH * (Y_NODE(8) - Y_MEAN) &
799 + AREA_TOP * (Z_NODE(8) - Z_MEAN) &
800 + AREA_WEST * (X_MEAN - X_NODE(1)) &
801 + AREA_SOUTH * (Y_MEAN - Y_NODE(1)) &
802 + AREA_BOTTOM * (Z_MEAN - Z_NODE(1)) &
803 + CUT_AREA * DEL_H
804
805 VOL_V(IJK) = VOLUME / 3.0D0
806 ENDIF
807
808 IF(IERROR==1) print*,'V_VEL IERROR',IJK
809 IF(IERROR==1) AXZ(IJK) = ZERO
810
811 IF(AXZ(IJK) <= TOL_SMALL_AREA*DX(I)*DZ(K)) THEN
812 WALL_V_AT(IJK) = .TRUE.
813 FLAG_N(IJK) = 0
814 NUMBER_OF_V_WALL_CELLS = NUMBER_OF_V_WALL_CELLS + 1
815 X_V(IJK) = CENTROID_CUT(1)
816 Y_V(IJK) = CENTROID_CUT(2)
817 Z_V(IJK) = CENTROID_CUT(3)
818 ENDIF
819
820 IF(N_QUADRIC>0) THEN
821 F_CUT_FACE = UNDEFINED
822 N_BC = 0
823 BC_V_ID(IJK) = 0
824 DO Q_ID = 1, N_QUADRIC
825 DO NODE = 1,N_CUT_FACE_NODES
826
827 CALL GET_F_QUADRIC(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
828 Q_ID,F_CUT_FACE(Q_ID),CLIP_FLAG)
829
830 IF(DABS(F_CUT_FACE(Q_ID)) < TOL_F) THEN
831 IF(BC_V_ID(IJK)/=BC_ID_Q(Q_ID)) N_BC = N_BC + 1
832 BC_V_ID(IJK) = BC_ID_Q(Q_ID)
833 ENDIF
834 ENDDO
835 IF(BC_V_ID(IJK)>0) THEN
836 IF(BC_TYPE_ENUM(BC_V_ID(IJK)) == CG_MI) EXIT
837
838 ENDIF
839
840 ENDDO
841
842 ENDIF
843
844 IF(N_POLYGON>0) THEN
845 DO NODE = 1,N_CUT_FACE_NODES
846 CALL EVAL_POLY_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
847 Q_ID,F_CUT_FACE(1),CLIP_FLAG,BCID)
848 IF(F_CUT_FACE(1)==ZERO) THEN
849 BC_V_ID(IJK) = BCID
850 IF(BC_V_ID(IJK)>0) THEN
851 IF(BC_TYPE_ENUM(BC_V_ID(IJK)) == CG_MI) EXIT
852 ENDIF
853 ENDIF
854 ENDDO
855 ENDIF
856
857
858 IF(N_USR_DEF>0) THEN
859 DO NODE = 1,N_CUT_FACE_NODES
860 CALL EVAL_USR_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
861 BCID,F_CUT_FACE(1),CLIP_FLAG)
862 IF(DABS(F_CUT_FACE(1)) < TOL_F) THEN
863 BC_V_ID(IJK) = BCID
864 IF(BC_V_ID(IJK)>0) THEN
865 IF(BC_TYPE_ENUM(BC_V_ID(IJK)) == CG_MI) EXIT
866 ENDIF
867 ENDIF
868 ENDDO
869 ENDIF
870
871
872 IF(USE_STL.OR.USE_MSH) THEN
873 DO NODE = 1,N_CUT_FACE_NODES
874 DO NN=1,N_FACET_AT(IJK)
875 NF=LIST_FACET_AT(IJK,NN)
876 CALL IS_POINT_INSIDE_FACET(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),&
877 COORD_CUT_FACE_NODES(3,NODE),NF,INSIDE_FACET)
878 IF(INSIDE_FACET) THEN
879 BC_V_ID(IJK) = BC_ID_STL_FACE(NF)
880 IF(BC_V_ID(IJK)>0) THEN
881 IF(BC_TYPE_ENUM(BC_V_ID(IJK)) == CG_MI) EXIT
882 ENDIF
883 ENDIF
884 ENDDO
885 ENDDO
886 ENDIF
887
888 CASE('W_MOMENTUM')
889 IF(I>ISTART1) AREA_WEST = AYZ_W(IMJK)
890 IF(J>JSTART1) AREA_SOUTH = AXZ_W(IJMK)
891 IF(K>KSTART1) AREA_BOTTOM = AXY_W(IJKM)
892
893 AYZ_W(IJK) = AREA_EAST
894 AXZ_W(IJK) = AREA_NORTH
895 AXY_W(IJK) = AREA_TOP
896 Area_W_CUT(IJK) = CUT_AREA
897
898 IF(AREA_EAST > ZERO) THEN
899 X_W_ec(IJK) = CENTROID_EAST(1)
900 Y_W_ec(IJK) = CENTROID_EAST(2)
901 Z_W_ec(IJK) = CENTROID_EAST(3)
902 ENDIF
903
904 IF(AREA_NORTH > ZERO) THEN
905 X_W_nc(IJK) = CENTROID_NORTH(1)
906 Y_W_nc(IJK) = CENTROID_NORTH(2)
907 Z_W_nc(IJK) = CENTROID_NORTH(3)
908 ENDIF
909
910 IF(AREA_TOP > ZERO) THEN
911 X_W_tc(IJK) = CENTROID_TOP(1)
912 Y_W_tc(IJK) = CENTROID_TOP(2)
913 Z_W_tc(IJK) = CENTROID_TOP(3)
914 ENDIF
915
916
917 (IJK,1) = AREA_EAST - AREA_WEST
918 NORMAL_W(IJK,2) = AREA_NORTH - AREA_SOUTH
919 NORMAL_W(IJK,3) = AREA_TOP - AREA_BOTTOM
920
921
922 NORM = sqrt(NORMAL_W(IJK,1)**2 + NORMAL_W(IJK,2)**2 + NORMAL_W(IJK,3)**2)
923
924 IF(NORM==ZERO) NORM = UNDEFINED
925
926 NORMAL_W(IJK,:) = NORMAL_W(IJK,:) / NORM
927
928 REFP_W(IJK,:) = CENTROID_CUT(:)
929
930
931 CALL GET_DEL_H(IJK,TYPE_OF_CELL,X_MEAN,Y_MEAN,Z_MEAN,Del_H,Nx,Ny,Nz)
932 IF(DEL_H == UNDEFINED) DEL_H = ZERO
933
934
935
936
937 = AREA_EAST * (X_NODE(8) - X_MEAN) &
938 + AREA_NORTH * (Y_NODE(8) - Y_MEAN) &
939 + AREA_TOP * (Z_NODE(8) - Z_MEAN) &
940 + AREA_WEST * (X_MEAN - X_NODE(1)) &
941 + AREA_SOUTH * (Y_MEAN - Y_NODE(1)) &
942 + AREA_BOTTOM * (Z_MEAN - Z_NODE(1)) &
943 + CUT_AREA * DEL_H
944
945 VOL_W(IJK) = VOLUME / 3.0D0
946
947 IF(IERROR==1) print*,'W_VEL IERROR',IJK
948 IF(IERROR==1) AXY(IJK) = ZERO
949
950
951 IF(AXY(IJK) <= TOL_SMALL_AREA*DX(I)*DY(J)) THEN
952 WALL_W_AT(IJK) = .TRUE.
953 FLAG_T(IJK) = 0
954 NUMBER_OF_W_WALL_CELLS = NUMBER_OF_W_WALL_CELLS + 1
955 X_W(IJK) = CENTROID_CUT(1)
956 Y_W(IJK) = CENTROID_CUT(2)
957 Z_W(IJK) = CENTROID_CUT(3)
958 ENDIF
959
960 IF(N_QUADRIC>0) THEN
961 F_CUT_FACE = UNDEFINED
962 N_BC = 0
963 BC_W_ID(IJK) = 0
964 DO Q_ID = 1, N_QUADRIC
965 DO NODE = 1,N_CUT_FACE_NODES
966 CALL GET_F_QUADRIC(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
967 Q_ID,F_CUT_FACE(Q_ID),CLIP_FLAG)
968 IF(DABS(F_CUT_FACE(Q_ID)) < TOL_F) THEN
969 IF(BC_W_ID(IJK)/=BC_ID_Q(Q_ID)) N_BC = N_BC + 1
970 BC_W_ID(IJK) = BC_ID_Q(Q_ID)
971 ENDIF
972 ENDDO
973 IF(BC_W_ID(IJK)>0) THEN
974 IF(BC_TYPE_ENUM(BC_W_ID(IJK)) == CG_MI) EXIT
975
976 ENDIF
977 ENDDO
978 ENDIF
979
980
981 IF(N_POLYGON>0) THEN
982 DO NODE = 1,N_CUT_FACE_NODES
983 CALL EVAL_POLY_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
984 Q_ID,F_CUT_FACE(1),CLIP_FLAG,BCID)
985 IF(F_CUT_FACE(1)==ZERO) THEN
986 BC_W_ID(IJK) = BCID
987 IF(BC_W_ID(IJK)>0) THEN
988 IF(BC_TYPE_ENUM(BC_W_ID(IJK)) == CG_MI) EXIT
989 ENDIF
990 ENDIF
991 ENDDO
992 ENDIF
993
994
995
996 IF(N_USR_DEF>0) THEN
997 DO NODE = 1,N_CUT_FACE_NODES
998 CALL EVAL_USR_FCT(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),COORD_CUT_FACE_NODES(3,NODE),&
999 BCID,F_CUT_FACE(1),CLIP_FLAG)
1000 IF(DABS(F_CUT_FACE(1)) < TOL_F) THEN
1001 BC_W_ID(IJK) = BCID
1002 IF(BC_W_ID(IJK)>0) THEN
1003 IF(BC_TYPE_ENUM(BC_W_ID(IJK)) == CG_MI) EXIT
1004 ENDIF
1005 ENDIF
1006 ENDDO
1007 ENDIF
1008
1009
1010 IF(USE_STL.OR.USE_MSH) THEN
1011 DO NODE = 1,N_CUT_FACE_NODES
1012 DO NN=1,N_FACET_AT(IJK)
1013 NF=LIST_FACET_AT(IJK,NN)
1014 CALL IS_POINT_INSIDE_FACET(COORD_CUT_FACE_NODES(1,NODE),COORD_CUT_FACE_NODES(2,NODE),&
1015 COORD_CUT_FACE_NODES(3,NODE),NF,INSIDE_FACET)
1016 IF(INSIDE_FACET) THEN
1017 BC_W_ID(IJK) = BC_ID_STL_FACE(NF)
1018 IF(BC_W_ID(IJK)>0) THEN
1019 IF(BC_TYPE_ENUM(BC_W_ID(IJK)) == CG_MI) EXIT
1020 ENDIF
1021 ENDIF
1022 ENDDO
1023 ENDDO
1024 ENDIF
1025
1026
1027 CASE DEFAULT
1028 WRITE(*,*)'SUBROUTINE: GET_CUT_CELL_VOLUME_AND_AREAS'
1029 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1030 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
1031 WRITE(*,*)'SCALAR'
1032 WRITE(*,*)'U_MOMENTUM'
1033 WRITE(*,*)'V_MOMENTUM'
1034 WRITE(*,*)'W_MOMENTUM'
1035 CALL MFIX_EXIT(myPE)
1036 END SELECT
1037
1038 RETURN
1039
1040 END SUBROUTINE GET_CUT_CELL_VOLUME_AND_AREAS
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055 SUBROUTINE GET_POLYGON_AREA_AND_CENTROID(NP,COORD,AREA,CENTROID,IERROR)
1056
1057 USE compar, only: mype
1058 USE cutcell
1059 USE exit, only: mfix_exit
1060 USE geometry, ONLY: no_k, zlength
1061 USE quadric, ONLY: cross_product
1062 USE param1, only: half, undefined, zero
1063
1064 IMPLICIT NONE
1065
1066 INTEGER, INTENT(INOUT) :: NP,IERROR
1067 DOUBLE PRECISION, INTENT(INOUT), DIMENSION(3) :: CENTROID
1068 DOUBLE PRECISION, INTENT(INOUT) :: AREA
1069 DOUBLE PRECISION, INTENT(INOUT), DIMENSION(3,15) :: COORD
1070
1071 INTEGER :: NN,NC,NEW_NP
1072 INTEGER :: LAST,LASTM
1073 DOUBLE PRECISION, DIMENSION(3,15) :: COORD_BCK
1074 DOUBLE PRECISION, DIMENSION(3) :: NORMAL,CP,SUMCP
1075 DOUBLE PRECISION, DIMENSION(3) :: LASTM_VEC,LAST_VEC
1076 DOUBLE PRECISION, DIMENSION(3) :: R,SUM_AR,AV
1077 DOUBLE PRECISION :: A,SUM_A
1078 LOGICAL, DIMENSION(15) :: KEEP
1079 DOUBLE PRECISION, DIMENSION(15) :: D
1080
1081
1082
1083
1084 IF(.FALSE.) THEN
1085 DO NN=1,NP
1086 D(NN) = sqrt(dot_product(COORD(:,NN),COORD(:,NN)))
1087 KEEP(NN) = .TRUE.
1088 COORD_BCK(:,NN) = COORD(:,NN)
1089 ENDDO
1090
1091 DO NN=1,NP-1
1092 DO NC=NN+1,NP
1093 IF(DABS(D(NN)-D(NC))<TOL_MERGE) KEEP(NC)=.FALSE.
1094 ENDDO
1095 ENDDO
1096
1097 NEW_NP = 0
1098
1099 DO NN=1,NP
1100 IF(KEEP(NN)) THEN
1101 NEW_NP = NEW_NP + 1
1102 COORD(:,NEW_NP) = COORD_BCK(:,NN)
1103 ENDIF
1104 ENDDO
1105
1106 NP = NEW_NP
1107 ENDIF
1108
1109 IF( NP < 2 ) THEN
1110 AREA = ZERO
1111 CENTROID = UNDEFINED
1112 RETURN
1113 ELSEIF( NP == 2 ) THEN
1114 IF(NO_K) THEN
1115 = sqrt((COORD(1,2)-COORD(1,1))**2 + (COORD(2,2)-COORD(2,1))**2) * ZLENGTH
1116 CENTROID(1) = HALF * (COORD(1,1)+COORD(1,2))
1117 CENTROID(2) = HALF * (COORD(2,1)+COORD(2,2))
1118 CENTROID(3) = ZERO
1119 RETURN
1120 ELSE
1121 AREA = ZERO
1122 CENTROID = UNDEFINED
1123 RETURN
1124
1125
1126
1127
1128
1129
1130 ENDIF
1131 ELSEIF( NP > 6 ) THEN
1132 WRITE(*,*)'CRITICAL ERROR: POLYGON WITH MORE THAN 6 POINTS.',NP
1133 DO NN=1,NP
1134 print*,COORD(:,NN)
1135 ENDDO
1136 WRITE(*,*)'MFiX will exit now.'
1137 CALL MFIX_EXIT(myPE)
1138 ENDIF
1139
1140
1141
1142
1143
1144
1145
1146 CALL REORDER_POLYGON(NP,COORD,NORMAL,IERROR)
1147
1148
1149
1150
1151
1152 (:,NP+1) = COORD(:,1)
1153
1154
1155
1156
1157
1158 = ZERO
1159
1160 DO NN = 1,NP
1161 CP = CROSS_PRODUCT(COORD(:,NN),COORD(:,NN+1))
1162 SUMCP = SUMCP + CP
1163 ENDDO
1164
1165
1166
1167
1168
1169 = DABS(HALF * DOT_PRODUCT(NORMAL, SUMCP))
1170
1171
1172
1173
1174
1175 = ZERO
1176 SUM_A = ZERO
1177 DO NN = 1,NP-2
1178 LAST = NN + 2
1179 LASTM = LAST - 1
1180 R = (COORD(:,1)+COORD(:,LASTM)+COORD(:,LAST))/3.0D0
1181 LASTM_VEC = COORD(:,LASTM)-COORD(:,1)
1182 LAST_VEC = COORD(:,LAST)-COORD(:,1)
1183 AV = CROSS_PRODUCT(LASTM_VEC,LAST_VEC)
1184 A = sqrt(dot_product(av(:),av(:)))
1185 SUM_AR = SUM_AR + A * R
1186 SUM_A = SUM_A + A
1187 ENDDO
1188 CENTROID = SUM_AR / SUM_A
1189
1190 RETURN
1191
1192 END SUBROUTINE GET_POLYGON_AREA_AND_CENTROID
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207 SUBROUTINE REORDER_POLYGON(NP,COORD,NORMAL,IERROR)
1208
1209 USE compar, ONLY: mype
1210 USE cutcell
1211 USE quadric, only: cross_product
1212 USE param1, only: zero
1213
1214 IMPLICIT NONE
1215 INTEGER :: I,K,NN,NP,NC,NEW_NP
1216 INTEGER :: LONGEST
1217 INTEGER :: I_SWAP,IERROR
1218
1219 DOUBLE PRECISION :: Xc,Yc,Zc,NORM,ANGLE,ANGLE_REF
1220 DOUBLE PRECISION, DIMENSION(3,15) :: COORD,COORD_BCK
1221 DOUBLE PRECISION, DIMENSION(3) :: NORMAL,SWAP
1222 DOUBLE PRECISION, DIMENSION(3) :: VECTMP,VECTMP2
1223
1224 LOGICAL, DIMENSION(15) :: KEEP
1225 DOUBLE PRECISION, DIMENSION(15) :: D
1226
1227
1228
1229
1230 IF(.FALSE.) THEN
1231 DO NN=1,NP
1232 D(NN) = sqrt(dot_product(COORD(:,NN),COORD(:,NN)))
1233 KEEP(NN) = .TRUE.
1234 COORD_BCK(:,NN) = COORD(:,NN)
1235 ENDDO
1236
1237 DO NN=1,NP-1
1238 DO NC=NN+1,NP
1239 IF(DABS(D(NN)-D(NC))<TOL_MERGE) KEEP(NC)=.FALSE.
1240 ENDDO
1241 ENDDO
1242
1243 NEW_NP = 0
1244
1245 DO NN=1,NP
1246 IF(KEEP(NN)) THEN
1247 NEW_NP = NEW_NP + 1
1248 COORD(:,NEW_NP) = COORD_BCK(:,NN)
1249 ENDIF
1250 ENDDO
1251
1252 NP = NEW_NP
1253 ENDIF
1254
1255
1256
1257
1258
1259 IF(NP<=2) RETURN
1260
1261
1262
1263
1264
1265 DO NN=1,NP
1266 ORDER(NN) = NN
1267 ENDDO
1268
1269
1270
1271
1272
1273 = ZERO
1274 Yc = ZERO
1275 Zc = ZERO
1276
1277 DO NN=1,NP
1278 Xc = Xc + COORD(1,NN)
1279 Yc = Yc + COORD(2,NN)
1280 Zc = Zc + COORD(3,NN)
1281 ENDDO
1282
1283 Xc = Xc / DBLE(NP)
1284 Yc = Yc / DBLE(NP)
1285 Zc = Zc / DBLE(NP)
1286
1287
1288
1289
1290
1291 = COORD(:,2)-COORD(:,1)
1292 VECTMP2 = COORD(:,3)-COORD(:,1)
1293 NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
1294
1295 NORM = sqrt(dot_product(NORMAL(:),NORMAL(:)))
1296 IF(NORM==ZERO) THEN
1297
1298
1299
1300
1301 = 1
1302
1303 RETURN
1304
1305 WRITE(*,*)'ERROR IN REORDER_POLYGON: NORMAL UNIT VECTOR HAS ZERO NORM'
1306 WRITE(*,*)'THIS USUALLY HAPPENS WHEN NODES WERE IMPROPERLY MERGED.'
1307 WRITE(*,*)'TO PREVENT THIS, INCREASE TOL_SNAP, OR DECREASE TOL_MERGE.'
1308 WRITE(*,*)'CURRENT VALUE OF TOL_SNAP = ', TOL_SNAP
1309 WRITE(*,*)'CURRENT VALUE OF TOL_MERGE = ', TOL_MERGE
1310 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1311 CALL MFIX_EXIT(MYPE)
1312 ELSE
1313 NORMAL = NORMAL / NORM
1314 ENDIF
1315
1316
1317
1318 = 3
1319 IF(ABS(NORMAL(1)) > ABS(NORMAL(2))) THEN
1320 IF(ABS(NORMAL(1)) > ABS(NORMAL(3))) LONGEST = 1
1321 ELSE
1322 IF(ABS(NORMAL(2)) > ABS(NORMAL(3))) LONGEST = 2
1323 ENDIF
1324
1325
1326
1327
1328 SELECT CASE(LONGEST)
1329
1330 CASE(1)
1331 DO I = 1, NP-1
1332 ANGLE_REF = ATAN2 (COORD(3,I) - Zc,COORD(2,I) - Yc)
1333 DO K = I+1,NP
1334 ANGLE = ATAN2 (COORD(3,K) - Zc,COORD(2,K) - Yc)
1335 IF(ANGLE < ANGLE_REF) THEN
1336 ANGLE_REF = ANGLE
1337 SWAP = COORD(:,K)
1338 COORD(:,K) = COORD(:,I)
1339 COORD(:,I) = SWAP
1340 I_SWAP = ORDER(K)
1341 ORDER(K) = ORDER(I)
1342 ORDER(I) = I_SWAP
1343 ENDIF
1344 END DO
1345 END DO
1346
1347 CASE(2)
1348 DO I = 1, NP-1
1349 ANGLE_REF = ATAN2 (COORD(1,I) - Xc,COORD(3,I) - Zc)
1350 DO K = I+1,NP
1351 ANGLE = ATAN2 (COORD(1,K) - Xc,COORD(3,K) - Zc)
1352 IF(ANGLE < ANGLE_REF) THEN
1353 ANGLE_REF = ANGLE
1354 SWAP = COORD(:,K)
1355 COORD(:,K) = COORD(:,I)
1356 COORD(:,I) = SWAP
1357 I_SWAP = ORDER(K)
1358 ORDER(K) = ORDER(I)
1359 ORDER(I) = I_SWAP
1360 ENDIF
1361 END DO
1362 END DO
1363
1364 CASE(3)
1365 DO I = 1, NP-1
1366 ANGLE_REF = ATAN2 (COORD(2,I) - Yc,COORD(1,I) - Xc)
1367 DO K = I+1,NP
1368 ANGLE = ATAN2 (COORD(2,K) - Yc,COORD(1,K) - Xc)
1369 IF(ANGLE < ANGLE_REF) THEN
1370 ANGLE_REF = ANGLE
1371 SWAP = COORD(:,K)
1372 COORD(:,K) = COORD(:,I)
1373 COORD(:,I) = SWAP
1374 I_SWAP = ORDER(K)
1375 ORDER(K) = ORDER(I)
1376 ORDER(I) = I_SWAP
1377 ENDIF
1378 END DO
1379 END DO
1380
1381 END SELECT
1382
1383 RETURN
1384
1385 END SUBROUTINE REORDER_POLYGON
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401 SUBROUTINE SORT(Array, NUM_ROWS,COL1,COL2)
1402
1403
1404 IMPLICIT NONE
1405
1406 INTEGER :: M,N,NUM_ROWS,COL,COL1,COL2,Saved_index
1407 DOUBLE PRECISION :: Swap,Local_min
1408 INTEGER, PARAMETER :: ARRAY_SIZE = 10000
1409 DOUBLE PRECISION, DIMENSION(ARRAY_SIZE,10) :: Array
1410
1411
1412 DO N = 1, NUM_ROWS-1
1413 Local_min = Array (N,COL1)
1414 Saved_index = N
1415 DO M = N,NUM_ROWS
1416 if (Array(M,COL1) < Local_min ) Then
1417 Local_min = Array (M,COL1)
1418 Saved_index = M
1419 endif
1420 END DO
1421
1422 DO COL = COL1,COL2
1423 Swap = Array(N,COL)
1424 Array(N,COL) = Array(Saved_index,COL)
1425 Array(Saved_index,COL) = Swap
1426 END DO
1427
1428 END DO
1429
1430 RETURN
1431
1432 END SUBROUTINE SORT
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447 SUBROUTINE GET_INTERPOLATION_TERMS_G(IJK,TYPE_OF_CELL,ALPHA_CUT,AW,HW,VELW)
1448
1449 USE bc
1450 USE compar, ONLY: mype
1451 USE cutcell
1452 USE exit, only: mfix_exit
1453 USE param1, only: one, zero, undefined
1454
1455 IMPLICIT NONE
1456 CHARACTER (LEN=*) :: TYPE_OF_CELL
1457 INTEGER :: IJK
1458 INTEGER :: BCV
1459 INTEGER :: BCT
1460 DOUBLE PRECISION :: ALPHA_CUT,AW,HW,VELW
1461
1462
1463
1464
1465
1466
1467
1468
1469 = ALPHA_CUT
1470 HW = ZERO
1471 VELW = ZERO
1472
1473 SELECT CASE (TYPE_OF_CELL)
1474
1475 CASE('U_MOMENTUM')
1476
1477 BCV = BC_U_ID(IJK)
1478 IF(BCV>0) THEN
1479 BCT = BC_TYPE_ENUM(BCV)
1480 SELECT CASE(BCT)
1481
1482 CASE(CG_NSW)
1483 AW = ALPHA_CUT
1484 HW = ZERO
1485 VELW = ZERO
1486
1487 CASE(CG_FSW)
1488 AW = ONE
1489 HW = ZERO
1490 VELW = ZERO
1491
1492 CASE(CG_PSW)
1493
1494 IF(BC_HW_G(BCV)==UNDEFINED) THEN
1495 = ALPHA_CUT
1496 HW = ZERO
1497 VELW = ZERO
1498 ELSEIF(BC_HW_G(BCV)==ZERO) THEN
1499 = ONE
1500 HW = ZERO
1501 VELW = ZERO
1502 ELSE
1503 = ONE
1504 = BC_HW_G(BCV)
1505 VELW = BC_UW_G(BCV)
1506 ENDIF
1507
1508 END SELECT
1509
1510 ENDIF
1511
1512
1513 CASE('V_MOMENTUM')
1514
1515 BCV = BC_V_ID(IJK)
1516 IF(BCV>0) THEN
1517 BCT = BC_TYPE_ENUM(BCV)
1518 SELECT CASE(BCT)
1519
1520 CASE(CG_NSW)
1521 AW = ALPHA_CUT
1522 HW = ZERO
1523 VELW = ZERO
1524
1525 CASE(CG_FSW)
1526 AW = ONE
1527 HW = ZERO
1528 VELW = ZERO
1529
1530 CASE(CG_PSW)
1531
1532 IF(BC_HW_G(BCV)==UNDEFINED) THEN
1533 = ALPHA_CUT
1534 HW = ZERO
1535 VELW = ZERO
1536 ELSEIF(BC_HW_G(BCV)==ZERO) THEN
1537 = ONE
1538 HW = ZERO
1539 VELW = ZERO
1540 ELSE
1541 = ONE
1542 = BC_HW_G(BCV)
1543 VELW = BC_VW_G(BCV)
1544 ENDIF
1545
1546 END SELECT
1547
1548 ENDIF
1549
1550
1551 CASE('W_MOMENTUM')
1552
1553 BCV = BC_W_ID(IJK)
1554 IF(BCV>0) THEN
1555 BCT = BC_TYPE_ENUM(BCV)
1556 SELECT CASE(BCT)
1557
1558 CASE(CG_NSW)
1559 AW = ALPHA_CUT
1560 HW = ZERO
1561 VELW = ZERO
1562
1563 CASE(CG_FSW)
1564 AW = ONE
1565 HW = ZERO
1566 VELW = ZERO
1567
1568 CASE(CG_PSW)
1569 IF(BC_HW_G(BCV)==UNDEFINED) THEN
1570 = ALPHA_CUT
1571 HW = ZERO
1572 VELW = ZERO
1573 ELSEIF(BC_HW_G(BCV)==ZERO) THEN
1574 = ONE
1575 HW = ZERO
1576 VELW = ZERO
1577 ELSE
1578 = ONE
1579 = BC_HW_G(BCV)
1580 VELW = BC_WW_G(BCV)
1581 ENDIF
1582
1583 END SELECT
1584
1585 ENDIF
1586
1587 CASE DEFAULT
1588 WRITE(*,*)'SUBROUTINE: GET_INTERPOLATION_TERMS_G'
1589 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1590 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
1591 WRITE(*,*)'U_MOMENTUM'
1592 WRITE(*,*)'V_MOMENTUM'
1593 WRITE(*,*)'W_MOMENTUM'
1594 CALL MFIX_EXIT(myPE)
1595 END SELECT
1596
1597 RETURN
1598
1599 END SUBROUTINE GET_INTERPOLATION_TERMS_G
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615 SUBROUTINE GET_INTERPOLATION_TERMS_S(IJK,M,TYPE_OF_CELL,ALPHA_CUT,AW,HW,VELW)
1616
1617 USE bc
1618 USE compar, ONLY: mype
1619 USE cutcell
1620 USE exit, only: mfix_exit
1621 USE param1, ONLY: one, zero, undefined
1622
1623 IMPLICIT NONE
1624 CHARACTER (LEN=*) :: TYPE_OF_CELL
1625 INTEGER :: IJK,M
1626 INTEGER :: BCV
1627 INTEGER :: BCT
1628 DOUBLE PRECISION :: ALPHA_CUT,AW,HW,VELW
1629
1630
1631
1632
1633
1634
1635
1636
1637 = ALPHA_CUT
1638 HW = ZERO
1639 VELW = ZERO
1640
1641 SELECT CASE (TYPE_OF_CELL)
1642
1643 CASE('U_MOMENTUM')
1644
1645 BCV = BC_U_ID(IJK)
1646 IF(BCV>0) THEN
1647 BCT = BC_TYPE_ENUM(BCV)
1648 SELECT CASE(BCT)
1649
1650 CASE(CG_NSW)
1651 AW = ALPHA_CUT
1652 HW = ZERO
1653 VELW = ZERO
1654
1655 CASE(CG_FSW)
1656 AW = ONE
1657 HW = ZERO
1658 VELW = ZERO
1659
1660 CASE(CG_PSW)
1661
1662 IF(BC_HW_S(BCV,M)==UNDEFINED) THEN
1663 = ALPHA_CUT
1664 HW = ZERO
1665 VELW = ZERO
1666 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN
1667 = ONE
1668 HW = ZERO
1669 VELW = ZERO
1670 ELSE
1671 = ALPHA_CUT
1672 HW = BC_HW_S(BCV,M)
1673 VELW = BC_UW_S(BCV,M)
1674 ENDIF
1675
1676 END SELECT
1677
1678 ENDIF
1679
1680
1681 CASE('V_MOMENTUM')
1682
1683 BCV = BC_V_ID(IJK)
1684 IF(BCV>0) THEN
1685 BCT = BC_TYPE_ENUM(BCV)
1686 SELECT CASE(BCT)
1687
1688 CASE(CG_NSW)
1689 AW = ALPHA_CUT
1690 HW = ZERO
1691 VELW = ZERO
1692
1693 CASE(CG_FSW)
1694 AW = ONE
1695 HW = ZERO
1696 VELW = ZERO
1697
1698 CASE(CG_PSW)
1699
1700 IF(BC_HW_S(BCV,M)==UNDEFINED) THEN
1701 = ALPHA_CUT
1702 HW = ZERO
1703 VELW = ZERO
1704 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN
1705 = ONE
1706 HW = ZERO
1707 VELW = ZERO
1708 ELSE
1709 = ALPHA_CUT
1710 HW = BC_HW_S(BCV,M)
1711 VELW = BC_VW_S(BCV,M)
1712 ENDIF
1713
1714 END SELECT
1715
1716 ENDIF
1717
1718
1719 CASE('W_MOMENTUM')
1720
1721 BCV = BC_W_ID(IJK)
1722 IF(BCV>0) THEN
1723 BCT = BC_TYPE_ENUM(BCV)
1724 SELECT CASE(BCT)
1725
1726 CASE(CG_NSW)
1727 AW = ALPHA_CUT
1728 HW = ZERO
1729 VELW = ZERO
1730
1731 CASE(CG_FSW)
1732 AW = ONE
1733 HW = ZERO
1734 VELW = ZERO
1735
1736 CASE(CG_PSW)
1737 IF(BC_HW_S(BCV,M)==UNDEFINED) THEN
1738 = ALPHA_CUT
1739 HW = ZERO
1740 VELW = ZERO
1741 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN
1742 = ONE
1743 HW = ZERO
1744 VELW = ZERO
1745 ELSE
1746 = ALPHA_CUT
1747 HW = BC_HW_S(BCV,M)
1748 VELW = BC_WW_S(BCV,M)
1749 ENDIF
1750
1751 END SELECT
1752
1753 ENDIF
1754
1755 CASE DEFAULT
1756 WRITE(*,*)'SUBROUTINE: GET_INTERPOLATION_TERMS_S'
1757 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
1758 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
1759 WRITE(*,*)'U_MOMENTUM'
1760 WRITE(*,*)'V_MOMENTUM'
1761 WRITE(*,*)'W_MOMENTUM'
1762 CALL MFIX_EXIT(myPE)
1763 END SELECT
1764
1765 RETURN
1766
1767 END SUBROUTINE GET_INTERPOLATION_TERMS_S
1768