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