File: N:\mfix\model\cartesian_grid\get_connectivity.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: GET_CONNECTIVITY                                       C
4     !  Purpose: Set flags for saclar cut cells, based on intersection      C
5     !  of the grid with the quadric(s)                                     C
6     !                                                                      C
7     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number #                                  Date: ##-###-##  C
11     !  Author: #                                                           C
12     !  Purpose: #                                                          C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15       SUBROUTINE GET_CONNECTIVITY(IJK,TYPE_OF_CELL,N_NEW_POINTS,N_NODES,CONNECT,X_NP,Y_NP,Z_NP,TOTAL_NUMBER_OF_INTERSECTIONS,&
16                  X_intersect,Y_intersect,Z_intersect)
17     
18           USE compar, ONLY: ijkend3
19           USE cutcell
20           USE cut_cell_preproc, ONLY: eval_f
21           USE functions, ONLY: FUNIJK
22           USE geometry, ONLY: DO_K, NO_K
23           USE indices, ONLY: I_OF, J_OF, K_OF
24           USE polygon, ONLY: n_polygon
25           USE quadric, ONLY: tol_f
26           USE param, ONLY: DIMENSION_3
27     
28           IMPLICIT NONE
29           INTEGER :: I,J,K,IM,JM,KM
30           INTEGER :: IJK,IMJK,IJMK,IJKM,IMJMK,IMJKM,IJMKM,IMJMKM
31           LOGICAL :: CLIP_FLAG
32           LOGICAL, DIMENSION(8) :: CORNER_INTERSECTION
33           INTEGER :: TOTAL_NUMBER_OF_INTERSECTIONS,NUMBER_OF_EDGE_INTERSECTIONS
34           INTEGER :: NUMBER_OF_CORNER_INTERSECTIONS,MAX_CORNER_INTERSECTIONS
35           INTEGER :: N_NODES,N_NEW_POINTS
36           INTEGER :: NODE,N_N1,N_N2,Q_ID,BCID
37           INTEGER, DIMENSION(DIMENSION_3,15) :: CONNECT
38           DOUBLE PRECISION, DIMENSION(DIMENSION_MAX_CUT_CELL) :: X_NP,Y_NP,Z_NP
39           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: X_intersect,Y_intersect,Z_intersect
40           CHARACTER (LEN=*) :: TYPE_OF_CELL
41     
42     !======================================================================
43     !  Get coordinates of eight nodes
44     !======================================================================
45     
46           CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
47     
48           I = I_OF(IJK)
49           J = J_OF(IJK)
50           K = K_OF(IJK)
51     
52           IM = I - 1
53           JM = J - 1
54           KM = K - 1
55     
56           IMJK   = FUNIJK(IM,J,K)
57           IJMK   = FUNIJK(I,JM,K)
58           IJKM   = FUNIJK(I,J,KM)
59     
60           IMJMK  = FUNIJK(IM,JM,K)
61           IMJKM  = FUNIJK(IM,J,KM)
62           IJMKM  = FUNIJK(I,JM,KM)
63     
64           IMJMKM = FUNIJK(IM,JM,KM)
65     
66     
67     !======================================================================
68     !  Evaluate Quadric at all corners
69     !======================================================================
70     
71              N_NODES = 0
72     
73              CORNER_INTERSECTION = .FALSE.
74              NUMBER_OF_CORNER_INTERSECTIONS = 0
75              NUMBER_OF_EDGE_INTERSECTIONS = 0
76     
77              IF(NO_K) THEN
78                 N_N1 = 5
79                 N_N2 = 8
80              ELSE
81                 N_N1 = 1
82                 N_N2 = 8
83              ENDIF
84     
85              DO NODE = N_N1,N_N2
86     
87                 Q_ID = 1
88                 CALL EVAL_F('QUADRIC',X_NODE(NODE),Y_NODE(NODE),Z_NODE(NODE),Q_ID,F_NODE(NODE),CLIP_FLAG)
89     
90                 CALL EVAL_F('POLYGON',X_NODE(NODE),Y_NODE(NODE),Z_NODE(NODE),N_POLYGON,F_NODE(NODE),CLIP_FLAG)
91     
92                 CALL EVAL_F('USR_DEF',X_NODE(NODE),Y_NODE(NODE),Z_NODE(NODE),N_USR_DEF,F_NODE(NODE),CLIP_FLAG)
93     
94                 CALL EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE,F_NODE(NODE),CLIP_FLAG,BCID)
95     
96                 IF (ABS(F_NODE(NODE)) < TOL_F ) THEN
97                    CORNER_INTERSECTION(NODE) = .TRUE.
98                    NUMBER_OF_CORNER_INTERSECTIONS = NUMBER_OF_CORNER_INTERSECTIONS + 1
99                    N_NODES = N_NODES + 1
100                    CONNECT(IJK,N_NODES) = IJK_OF_NODE(NODE)
101                 ENDIF
102     
103     
104                 IF(SNAP(IJK_OF_NODE(NODE))) THEN
105                    CORNER_INTERSECTION(NODE) = .TRUE.
106                    NUMBER_OF_CORNER_INTERSECTIONS = NUMBER_OF_CORNER_INTERSECTIONS + 1
107                    N_NODES = N_NODES + 1
108                    CONNECT(IJK,N_NODES) = IJK_OF_NODE(NODE)
109                 ENDIF
110     
111              END DO
112     
113     !======================================================================
114     !  Count the number of edge intersections (excluding corner intersections)
115     !  For each new edge intersection found:
116     !     - Increment the total number of points
117     !     - Store the location of the additional point
118     !======================================================================
119     
120                 IF(DO_K) THEN
121     
122                    IF(INTERSECT_X(IJMKM)) THEN  ! Edge 1  = Nodes 1-2
123                       IF((.NOT.CORNER_INTERSECTION(1)).AND.(.NOT.CORNER_INTERSECTION(2))) THEN
124                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
125                          N_NEW_POINTS = N_NEW_POINTS + 1
126                          X_NP(N_NEW_POINTS) = X_intersect(IJMKM)
127                          Y_NP(N_NEW_POINTS) = Y_NODE(1)
128                          Z_NP(N_NEW_POINTS) = Z_NODE(1)
129                          N_NODES = N_NODES + 1
130                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
131                       ENDIF
132                    ENDIF
133     
134                    IF(INTERSECT_Y(IJKM)) THEN  ! Edge 2  = Nodes 2-4
135                       IF((.NOT.CORNER_INTERSECTION(2)).AND.(.NOT.CORNER_INTERSECTION(4))) THEN
136                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
137                          N_NEW_POINTS = N_NEW_POINTS + 1
138                          X_NP(N_NEW_POINTS) = X_NODE(2)
139                          Y_NP(N_NEW_POINTS) = Y_intersect(IJKM)
140                          Z_NP(N_NEW_POINTS) = Z_NODE(2)
141                          N_NODES = N_NODES + 1
142                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
143                       ENDIF
144                    ENDIF
145     
146                    IF(INTERSECT_X(IJKM)) THEN  ! Edge 3  = Nodes 3-4
147                       IF((.NOT.CORNER_INTERSECTION(3)).AND.(.NOT.CORNER_INTERSECTION(4))) THEN
148                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
149                          N_NEW_POINTS = N_NEW_POINTS + 1
150                          X_NP(N_NEW_POINTS) = X_intersect(IJKM)
151                          Y_NP(N_NEW_POINTS) = Y_NODE(3)
152                          Z_NP(N_NEW_POINTS) = Z_NODE(3)
153                          N_NODES = N_NODES + 1
154                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
155                       ENDIF
156                    ENDIF
157     
158                    IF(INTERSECT_Y(IMJKM)) THEN  ! Edge 4  = Nodes 1-3
159                       IF((.NOT.CORNER_INTERSECTION(1)).AND.(.NOT.CORNER_INTERSECTION(3))) THEN
160                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
161                          N_NEW_POINTS = N_NEW_POINTS + 1
162                          X_NP(N_NEW_POINTS) = X_NODE(1)
163                          Y_NP(N_NEW_POINTS) = Y_intersect(IMJKM)
164                          Z_NP(N_NEW_POINTS) = Z_NODE(1)
165                          N_NODES = N_NODES + 1
166                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
167                       ENDIF
168                    ENDIF
169     
170                 ENDIF
171     
172     
173                 IF(INTERSECT_X(IJMK)) THEN  ! Edge 5  = Nodes 5-6
174                    IF((.NOT.CORNER_INTERSECTION(5)).AND.(.NOT.CORNER_INTERSECTION(6))) THEN
175                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
176                       N_NEW_POINTS = N_NEW_POINTS + 1
177                       X_NP(N_NEW_POINTS) = X_intersect(IJMK)
178                       Y_NP(N_NEW_POINTS) = Y_NODE(5)
179                       Z_NP(N_NEW_POINTS) = Z_NODE(5)
180                       N_NODES = N_NODES + 1
181                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
182                    ENDIF
183                 ENDIF
184     
185                 IF(INTERSECT_Y(IJK)) THEN  ! Edge 6  = Nodes 6-8
186                    IF((.NOT.CORNER_INTERSECTION(6)).AND.(.NOT.CORNER_INTERSECTION(8))) THEN
187                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
188                       N_NEW_POINTS = N_NEW_POINTS + 1
189                       X_NP(N_NEW_POINTS) = X_NODE(6)
190                       Y_NP(N_NEW_POINTS) = Y_intersect(IJK)
191                       Z_NP(N_NEW_POINTS) = Z_NODE(6)
192                       N_NODES = N_NODES + 1
193                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
194                    ENDIF
195                 ENDIF
196     
197                 IF(INTERSECT_X(IJK)) THEN  ! Edge 7  = Nodes 7-8
198                    IF((.NOT.CORNER_INTERSECTION(7)).AND.(.NOT.CORNER_INTERSECTION(8))) THEN
199                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
200                       N_NEW_POINTS = N_NEW_POINTS + 1
201                       X_NP(N_NEW_POINTS) = X_intersect(IJK)
202                       Y_NP(N_NEW_POINTS) = Y_NODE(7)
203                       Z_NP(N_NEW_POINTS) = Z_NODE(7)
204                       N_NODES = N_NODES + 1
205                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
206                    ENDIF
207                 ENDIF
208     
209                 IF(INTERSECT_Y(IMJK)) THEN  ! Edge 8  = Nodes 5-7
210                    IF((.NOT.CORNER_INTERSECTION(5)).AND.(.NOT.CORNER_INTERSECTION(7))) THEN
211                       NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
212                       N_NEW_POINTS = N_NEW_POINTS + 1
213                       X_NP(N_NEW_POINTS) = X_NODE(5)
214                       Y_NP(N_NEW_POINTS) = Y_intersect(IMJK)
215                       Z_NP(N_NEW_POINTS) = Z_NODE(5)
216                       N_NODES = N_NODES + 1
217                       CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
218                    ENDIF
219                 ENDIF
220     
221     
222                 IF(DO_K) THEN
223     
224                    IF(INTERSECT_Z(IMJMK)) THEN  ! Edge 9  = Nodes 1-5
225                       IF((.NOT.CORNER_INTERSECTION(1)).AND.(.NOT.CORNER_INTERSECTION(5))) THEN
226                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
227                          N_NEW_POINTS = N_NEW_POINTS + 1
228                          X_NP(N_NEW_POINTS) = X_NODE(1)
229                          Y_NP(N_NEW_POINTS) = Y_NODE(1)
230                          Z_NP(N_NEW_POINTS) = Z_intersect(IMJMK)
231                          N_NODES = N_NODES + 1
232                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
233                       ENDIF
234                    ENDIF
235     
236                    IF(INTERSECT_Z(IJMK))  THEN  ! Edge 10 = Nodes 2-6
237                       IF((.NOT.CORNER_INTERSECTION(2)).AND.(.NOT.CORNER_INTERSECTION(6))) THEN
238                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
239                          N_NEW_POINTS = N_NEW_POINTS + 1
240                          X_NP(N_NEW_POINTS) = X_NODE(2)
241                          Y_NP(N_NEW_POINTS) = Y_NODE(2)
242                          Z_NP(N_NEW_POINTS) = Z_intersect(IJMK)
243                          N_NODES = N_NODES + 1
244                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
245                       ENDIF
246                    ENDIF
247     
248                    IF(INTERSECT_Z(IJK))  THEN  ! Edge 11 = Nodes 4-8
249                       IF((.NOT.CORNER_INTERSECTION(4)).AND.(.NOT.CORNER_INTERSECTION(8))) THEN
250                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
251                          N_NEW_POINTS = N_NEW_POINTS + 1
252                          X_NP(N_NEW_POINTS) = X_NODE(4)
253                          Y_NP(N_NEW_POINTS) = Y_NODE(4)
254                          Z_NP(N_NEW_POINTS) = Z_intersect(IJK)
255                          N_NODES = N_NODES + 1
256                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
257                       ENDIF
258                    ENDIF
259     
260                    IF(INTERSECT_Z(IMJK)) THEN  ! Edge 12 = Nodes 3-7
261                       IF((.NOT.CORNER_INTERSECTION(3)).AND.(.NOT.CORNER_INTERSECTION(7))) THEN
262                          NUMBER_OF_EDGE_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + 1
263                          N_NEW_POINTS = N_NEW_POINTS + 1
264                          X_NP(N_NEW_POINTS) = X_NODE(3)
265                          Y_NP(N_NEW_POINTS) = Y_NODE(3)
266                          Z_NP(N_NEW_POINTS) = Z_intersect(IMJK)
267                          N_NODES = N_NODES + 1
268                          CONNECT(IJK,N_NODES) = N_NEW_POINTS + IJKEND3
269                       ENDIF
270                    ENDIF
271     
272                 ENDIF
273     
274     !======================================================================
275     !  Count the total number of intersections (corner and edge intersections)
276     !======================================================================
277     
278                 IF(NO_K) THEN
279                    MAX_CORNER_INTERSECTIONS = 2
280                 ELSE
281                    MAX_CORNER_INTERSECTIONS = 4
282                 ENDIF
283     
284     
285                IF(NUMBER_OF_CORNER_INTERSECTIONS > MAX_CORNER_INTERSECTIONS) THEN
286                   IF(PRINT_WARNINGS) THEN
287                      WRITE(*,*)'WARNING:',NUMBER_OF_CORNER_INTERSECTIONS,&
288                             ' CORNER INTERSECTIONS DETECTED IN CELL IJK=',IJK
289                      WRITE(*,*)'THIS USUALLY INDICATE A FALSE CUT-CELL (CORNER CELL)'
290     !                 WRITE(*,*)'RESETTING NUMBER_OF_CORNER_INTERSECTIONS TO', MAX_CORNER_INTERSECTIONS
291                      WRITE(*,*)'REMOVING CUT CELL'
292                   ENDIF
293     !              NUMBER_OF_CORNER_INTERSECTIONS = MAX_CORNER_INTERSECTIONS
294     
295                   ! Force the total number of intersections to be zero, and therefore, the cell will be considered as a non-cut cell
296     !              NUMBER_OF_CORNER_INTERSECTIONS = -NUMBER_OF_EDGE_INTERSECTIONS
297     
298                   ! Force the total number of intersections to be -one, and therefore, the cell will be considered as a non-cut cell
299                   NUMBER_OF_CORNER_INTERSECTIONS = -NUMBER_OF_EDGE_INTERSECTIONS -1
300     
301     
302                ENDIF
303     
304                 TOTAL_NUMBER_OF_INTERSECTIONS = NUMBER_OF_EDGE_INTERSECTIONS + NUMBER_OF_CORNER_INTERSECTIONS
305     
306           RETURN
307     
308           END SUBROUTINE GET_CONNECTIVITY
309     
310     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
311     !                                                                      C
312     !  Module name: GET_CELL_NODE_COORDINATES                              C
313     !  Purpose: Get the cell corners (x,y,z) coordinates                   C
314     !                                                                      C
315     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
316     !  Reviewer:                                          Date:            C
317     !                                                                      C
318     !  Revision Number #                                  Date: ##-###-##  C
319     !  Author: #                                                           C
320     !  Purpose: #                                                          C
321     !                                                                      C
322     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
323       SUBROUTINE GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
324     
325           USE compar, ONLY: mype
326           USE cutcell
327           USE exit, only: mfix_exit
328           USE functions, ONLY: FUNIJK
329           USE geometry, ONLY: NO_K, dx, dy, dz
330           USE indices, ONLY: I_OF, J_OF, K_OF
331           USE param1, ONLY: HALF, ZERO
332     
333           IMPLICIT NONE
334           CHARACTER (LEN=*) :: TYPE_OF_CELL
335           DOUBLE PRECISION:: Xw,Xe,Yn,Ys,Zb,Zt
336           INTEGER :: I,J,K,IP,JP,KP,IM,JM,KM
337           INTEGER :: IJK,IMJK,IJMK,IJKM,IMJMK,IMJKM,IJMKM,IMJMKM
338     
339           I = I_OF(IJK)
340           J = J_OF(IJK)
341           K = K_OF(IJK)
342     
343           IP = I + 1
344           JP = J + 1
345           KP = K + 1
346     
347           IM = I - 1
348           JM = J - 1
349           KM = K - 1
350     
351           IMJK   = FUNIJK(IM,J,K)
352           IJMK   = FUNIJK(I,JM,K)
353           IJKM   = FUNIJK(I,J,KM)
354     
355           IMJMK  = FUNIJK(IM,JM,K)
356           IMJKM  = FUNIJK(IM,J,KM)
357           IJMKM  = FUNIJK(I,JM,KM)
358     
359           IMJMKM = FUNIJK(IM,JM,KM)
360     
361           IJK_OF_NODE(0) = IJK
362           IJK_OF_NODE(1) = IMJMKM
363           IJK_OF_NODE(2) = IJMKM
364           IJK_OF_NODE(3) = IMJKM
365           IJK_OF_NODE(4) = IJKM
366           IJK_OF_NODE(5) = IMJMK
367           IJK_OF_NODE(6) = IJMK
368           IJK_OF_NODE(7) = IMJK
369           IJK_OF_NODE(8) = IJK
370           IJK_OF_NODE(15) = IJK
371     
372     
373           SELECT CASE (TYPE_OF_CELL)
374              CASE('SCALAR')
375                 Xw = XG_E(I) - DX(I)          ! west face location
376                 Xe = XG_E(I)                  ! east face location
377     
378                 Ys = YG_N(J) - DY(J)          ! south face location
379                 Yn = YG_N(J)                  ! north face location
380     
381                 IF(NO_K) THEN
382                    Zb = ZERO                  ! bottom face location
383                    Zt = ZERO                  ! top face location
384                 ELSE
385                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
386                    Zt = ZG_T(K)               ! top face location
387                 ENDIF
388     
389              CASE('U_MOMENTUM')
390                 Xw = XG_E(I) - HALF * DX(I)   ! west face location
391                 Xe = XG_E(I) + HALF * DX(IP)  ! east face location
392     
393                 Ys = YG_N(J) - DY(J)          ! south face location
394                 Yn = YG_N(J)                  ! north face location
395     
396                 IF(NO_K) THEN
397                    Zb = ZERO                  ! bottom face location
398                    Zt = ZERO                  ! top face location
399                 ELSE
400                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
401                    Zt = ZG_T(K)               ! top face location
402                 ENDIF
403     
404              CASE('V_MOMENTUM')
405                 Xw = XG_E(I) - DX(I)          ! west face location
406                 Xe = XG_E(I)                  ! east face location
407     
408                 Ys = YG_N(J) - HALF * DY(J)   ! south face location
409                 Yn = YG_N(J) + HALF * DY(JP)  ! north face location
410     
411                 IF(NO_K) THEN
412                    Zb = ZERO                  ! bottom face location
413                    Zt = ZERO                  ! top face location
414                 ELSE
415                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
416                    Zt = ZG_T(K)               ! top face location
417                 ENDIF
418     
419              CASE('W_MOMENTUM')
420                 Xw = XG_E(I) - DX(I)          ! west face location
421                 Xe = XG_E(I)                  ! east face location
422     
423                 Ys = YG_N(J) - DY(J)          ! south face location
424                 Yn = YG_N(J)                  ! north face location
425     
426                 Zb = ZG_T(K) - HALF * DZ(K)   ! bottom face location
427                 Zt = ZG_T(K) + HALF * DZ(KP)  ! top face location
428     
429     
430              CASE DEFAULT
431                 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
432                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
433                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
434                 WRITE(*,*)'SCALAR'
435                 WRITE(*,*)'U_MOMENTUM'
436                 WRITE(*,*)'V_MOMENTUM'
437                 WRITE(*,*)'W_MOMENTUM'
438                 CALL MFIX_EXIT(myPE)
439           END SELECT
440     
441     !     CELL CENTER :
442           X_NODE(0) = HALF * ( xw + xe )
443           Y_NODE(0) = HALF * ( ys + yn )
444           Z_NODE(0) = HALF * ( zb + zt )
445     
446     !     NODE 1 : IM,JM,KM
447           X_NODE(1) = xw
448           Y_NODE(1) = ys
449           Z_NODE(1) = zb
450     
451     !     NODE 2 : I,JM,KM
452           X_NODE(2) = xe
453           Y_NODE(2) = ys
454           Z_NODE(2) = zb
455     
456     !     NODE 3 : IM,J,KM
457           X_NODE(3) = xw
458           Y_NODE(3) = yn
459           Z_NODE(3) = zb
460     
461     !     NODE 4 : I,J,KM
462           X_NODE(4) = xe
463           Y_NODE(4) = yn
464           Z_NODE(4) = zb
465     
466     !     NODE 5 : IM,JM,K
467           X_NODE(5) = xw
468           Y_NODE(5) = ys
469           Z_NODE(5) = zt
470     
471     !     NODE 6 : I,JM,K
472           X_NODE(6) = xe
473           Y_NODE(6) = ys
474           Z_NODE(6) = zt
475     
476     !     NODE 7 : IM,J,K
477           X_NODE(7) = xw
478           Y_NODE(7) = yn
479           Z_NODE(7) = zt
480     
481     !     NODE 8 : I,J,K
482           X_NODE(8) = xe
483           Y_NODE(8) = yn
484           Z_NODE(8) = zt
485     
486           RETURN
487     
488           END SUBROUTINE GET_CELL_NODE_COORDINATES
489     
490     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
491     !                                                                      C
492     !  Module name: GET_GLOBAL_CELL_NODE_COORDINATES                       C
493     !  Purpose: Get the cell corners (x,y,z) coordinates                   C
494     !                                                                      C
495     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
496     !  Reviewer:                                          Date:            C
497     !                                                                      C
498     !  Revision Number #                                  Date: ##-###-##  C
499     !  Author: #                                                           C
500     !  Purpose: #                                                          C
501     !                                                                      C
502     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
503       SUBROUTINE GET_GLOBAL_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
504     
505           USE compar, ONLY: MYPE
506           USE cutcell
507           USE exit, only: mfix_exit
508           USE functions, ONLY: FUNIJK_GL
509           USE geometry, ONLY: NO_K, dx, dy, dz
510           USE param1, only: half, zero
511           USE vtk, ONLY: GLOBAL_I_OF, GLOBAL_J_OF, GLOBAL_K_OF
512     
513           IMPLICIT NONE
514           CHARACTER (LEN=*) :: TYPE_OF_CELL
515           DOUBLE PRECISION:: Xw,Xe,Yn,Ys,Zb,Zt
516           INTEGER :: I,J,K,IP,JP,KP,IM,JM,KM
517           INTEGER :: IJK,IMJK,IJMK,IJKM,IMJMK,IMJKM,IJMKM,IMJMKM
518     
519           I = GLOBAL_I_OF(IJK)
520           J = GLOBAL_J_OF(IJK)
521           K = GLOBAL_K_OF(IJK)
522     
523           IP = I + 1
524           JP = J + 1
525           KP = K + 1
526     
527           IM = I - 1
528           JM = J - 1
529           KM = K - 1
530     
531           IMJK   = FUNIJK_GL(IM,J,K)
532           IJMK   = FUNIJK_GL(I,JM,K)
533           IJKM   = FUNIJK_GL(I,J,KM)
534     
535           IMJMK  = FUNIJK_GL(IM,JM,K)
536           IMJKM  = FUNIJK_GL(IM,J,KM)
537           IJMKM  = FUNIJK_GL(I,JM,KM)
538     
539           IMJMKM = FUNIJK_GL(IM,JM,KM)
540     
541           IJK_OF_NODE(1) = IMJMKM
542           IJK_OF_NODE(2) = IJMKM
543           IJK_OF_NODE(3) = IMJKM
544           IJK_OF_NODE(4) = IJKM
545           IJK_OF_NODE(5) = IMJMK
546           IJK_OF_NODE(6) = IJMK
547           IJK_OF_NODE(7) = IMJK
548           IJK_OF_NODE(8) = IJK
549     
550           SELECT CASE (TYPE_OF_CELL)
551              CASE('SCALAR')
552                 Xw = XG_E(I) - DX(I)          ! west face location
553                 Xe = XG_E(I)                  ! east face location
554     
555                 Ys = YG_N(J) - DY(J)          ! south face location
556                 Yn = YG_N(J)                  ! north face location
557     
558                 IF(NO_K) THEN
559                    Zb = ZERO                  ! bottom face location
560                    Zt = ZERO                  ! top face location
561                 ELSE
562                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
563                    Zt = ZG_T(K)               ! top face location
564                 ENDIF
565     
566              CASE('U_MOMENTUM')
567                 Xw = XG_E(I) - HALF * DX(I)   ! west face location
568                 Xe = XG_E(I) + HALF * DX(IP)  ! east face location
569     
570                 Ys = YG_N(J) - DY(J)          ! south face location
571                 Yn = YG_N(J)                  ! north face location
572     
573                 IF(NO_K) THEN
574                    Zb = ZERO                  ! bottom face location
575                    Zt = ZERO                  ! top face location
576                 ELSE
577                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
578                    Zt = ZG_T(K)               ! top face location
579                 ENDIF
580     
581              CASE('V_MOMENTUM')
582                 Xw = XG_E(I) - DX(I)          ! west face location
583                 Xe = XG_E(I)                  ! east face location
584     
585                 Ys = YG_N(J) - HALF * DY(J)   ! south face location
586                 Yn = YG_N(J) + HALF * DY(JP)  ! north face location
587     
588                 IF(NO_K) THEN
589                    Zb = ZERO                  ! bottom face location
590                    Zt = ZERO                  ! top face location
591                 ELSE
592                    Zb = ZG_T(K) - DZ(K)       ! bottom face location
593                    Zt = ZG_T(K)               ! top face location
594                 ENDIF
595     
596              CASE('W_MOMENTUM')
597                 Xw = XG_E(I) - DX(I)          ! west face location
598                 Xe = XG_E(I)                  ! east face location
599     
600                 Ys = YG_N(J) - DY(J)          ! south face location
601                 Yn = YG_N(J)                  ! north face location
602     
603                 Zb = ZG_T(K) - HALF * DZ(K)   ! bottom face location
604                 Zt = ZG_T(K) + HALF * DZ(KP)  ! top face location
605     
606     
607              CASE DEFAULT
608                 WRITE(*,*)'SUBROUTINE: GET_CELL_NODE_COORDINATES'
609                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
610                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
611                 WRITE(*,*)'SCALAR'
612                 WRITE(*,*)'U_MOMENTUM'
613                 WRITE(*,*)'V_MOMENTUM'
614                 WRITE(*,*)'W_MOMENTUM'
615                 CALL MFIX_EXIT(myPE)
616           END SELECT
617     
618     !     CELL CENTER :
619           X_NODE(0) = HALF * ( xw + xe )
620           Y_NODE(0) = HALF * ( ys + yn )
621           Z_NODE(0) = HALF * ( zb + zt )
622     
623     !     NODE 1 : IM,JM,KM
624           X_NODE(1) = xw
625           Y_NODE(1) = ys
626           Z_NODE(1) = zb
627     
628     !     NODE 2 : I,JM,KM
629           X_NODE(2) = xe
630           Y_NODE(2) = ys
631           Z_NODE(2) = zb
632     
633     !     NODE 3 : IM,J,KM
634           X_NODE(3) = xw
635           Y_NODE(3) = yn
636           Z_NODE(3) = zb
637     
638     !     NODE 4 : I,J,KM
639           X_NODE(4) = xe
640           Y_NODE(4) = yn
641           Z_NODE(4) = zb
642     
643     !     NODE 5 : IM,JM,K
644           X_NODE(5) = xw
645           Y_NODE(5) = ys
646           Z_NODE(5) = zt
647     
648     !     NODE 6 : I,JM,K
649           X_NODE(6) = xe
650           Y_NODE(6) = ys
651           Z_NODE(6) = zt
652     
653     !     NODE 7 : IM,J,K
654           X_NODE(7) = xw
655           Y_NODE(7) = yn
656           Z_NODE(7) = zt
657     
658     !     NODE 8 : I,J,K
659           X_NODE(8) = xe
660           Y_NODE(8) = yn
661           Z_NODE(8) = zt
662     
663           RETURN
664     
665           END SUBROUTINE GET_GLOBAL_CELL_NODE_COORDINATES
666