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