File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/get_poly_data.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 SUBROUTINE GET_POLY_DATA
26
27
28
29
30
31 USE param
32 USE param1
33 USE physprop
34 USE fldvar
35 USE run
36 USE scalars
37 USE funits
38 USE rxns
39 USE compar
40 USE mpi_utility
41 USE progress_bar
42 USE polygon
43 IMPLICIT NONE
44
45 INTEGER :: POLY,V,N,NSKIP
46 LOGICAL :: PRESENT
47
48
49 WRITE(*,2000) 'READING polygon geometry from poly.dat...'
50
51 INQUIRE(FILE='poly.dat',EXIST=PRESENT)
52 IF(.NOT.PRESENT) THEN
53 IF(MyPE == PE_IO) THEN
54 WRITE(*,"('(PE ',I3,'): input data file, ',A11,' is missing: run aborted')") &
55 myPE,'poly.dat'
56 ENDIF
57 CALL MFIX_EXIT(MYPE)
58 ENDIF
59
60
61
62
63 OPEN(UNIT=333, FILE='poly.dat', STATUS='OLD', ERR=910)
64
65 NSKIP = 13
66 DO N=1,NSKIP
67 READ(333,*,ERR=920,END=930)
68 ENDDO
69
70 READ(333,*,ERR=920,END=930) N_POLYGON
71
72 DO POLY = 1, N_POLYGON
73 READ(333,*,ERR=920,END=930) N_VERTEX(POLY), POLY_SIGN(POLY)
74 IF(DABS(POLY_SIGN(POLY))/=ONE) THEN
75 WRITE(*,*)'ERROR WHILE READIND poly.dat:'
76 WRITE(*,*)'POLYGON SIGN MUST BE +1.0 or -1.0.'
77 CALL MFIX_EXIT(myPE)
78 ENDIF
79 DO V = 1,N_VERTEX(POLY)
80 READ(333,*,ERR=920,END=930) X_VERTEX(POLY,V),Y_VERTEX(POLY,V),BC_ID_P(POLY,V)
81 ENDDO
82 X_VERTEX(POLY,N_VERTEX(POLY) + 1) = X_VERTEX(POLY, 1)
83 (POLY,N_VERTEX(POLY) + 1) = Y_VERTEX(POLY, 1)
84 ENDDO
85
86 CLOSE(333)
87
88 WRITE(*,2010) 'Polygon geometry successfully read for ', N_POLYGON, ' polygon(s).'
89
90 RETURN
91
92
93
94
95
96 CONTINUE
97 WRITE (*, 1500)
98 CALL MFIX_EXIT(myPE)
99 920 CONTINUE
100 WRITE (*, 1600)
101 CALL MFIX_EXIT(myPE)
102 930 CONTINUE
103 WRITE (*, 1700)
104 CALL MFIX_EXIT(myPE)
105
106 FORMAT(/1X,70('*')//' From: GET_POLY_DATA',/' Message: ',&
107 'Unable to open poly.dat file',/1X,70('*')/)
108 1600 FORMAT(/1X,70('*')//' From: GET_POLY_DATA',/' Message: ',&
109 'Error while reading poly.dat file',/1X,70('*')/)
110 1700 FORMAT(/1X,70('*')//' From: GET_POLY_DATA',/' Message: ',&
111 'End of file reached while reading poly.dat file',/1X,70('*')/)
112 2000 FORMAT(1X,A)
113 2010 FORMAT(1X,A,I4,A)
114
115 END SUBROUTINE GET_POLY_DATA
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131 SUBROUTINE EVAL_POLY_FCT(x1,x2,x3,Q,f_pol,CLIP_FLAG,BCID)
132
133 USE param
134 USE param1
135 USE parallel
136 USE constant
137 USE run
138 USE toleranc
139 USE geometry
140 USE indices
141 USE compar
142 USE sendrecv
143 USE fldvar
144 USE quadric
145 USE cutcell
146 USE polygon
147
148 IMPLICIT NONE
149 DOUBLE PRECISION x1,x2,x3
150 DOUBLE PRECISION f_pol
151 DOUBLE PRECISION, DIMENSION(DIM_POLYGON) :: F_POLY
152 INTEGER :: Q
153 LOGICAL :: CLIP_X,CLIP_Y,CLIP_FLAG
154 LOGICAL :: test1,test2
155
156 INTEGER :: POLY,N_EDGE,COUNTER,EDGE
157 DOUBLE PRECISION :: X_VERTEX_MIN,X_VERTEX_MAX,Y_VERTEX_MIN,Y_VERTEX_MAX
158 DOUBLE PRECISION :: XV1,XV2,YV1,YV2
159 DOUBLE PRECISION :: x_west,x_east,y_north,y_south,slope,x_star,D1,D2
160 DOUBLE PRECISION :: max_slope = 200.0D0
161 DOUBLE PRECISION :: TOL_xstar
162 INTEGER :: BCID
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184 IF(N_POLYGON < 1) RETURN
185
186 BCID = -1
187
188 DO POLY = 1,N_POLYGON
189
190 N_EDGE = N_VERTEX(POLY)
191
192 X_VERTEX_MIN = MINVAL(X_VERTEX(POLY,:)) - 2.0D0 * TOL_POLY
193 X_VERTEX_MAX = MAXVAL(X_VERTEX(POLY,:)) + 2.0D0 * TOL_POLY
194
195 Y_VERTEX_MIN = MINVAL(Y_VERTEX(POLY,:)) - 2.0D0 * TOL_POLY
196 Y_VERTEX_MAX = MAXVAL(Y_VERTEX(POLY,:)) + 2.0D0 * TOL_POLY
197
198 CLIP_X = (X_VERTEX_MIN <= x1).AND.( x1 <= X_VERTEX_MAX)
199 CLIP_Y = (Y_VERTEX_MIN <= x2).AND.( x2 <= Y_VERTEX_MAX)
200
201
202 IF(CLIP_X.AND.CLIP_Y) THEN
203
204 COUNTER = 0
205
206 DO EDGE = 1, N_EDGE
207
208 XV1 = X_VERTEX(POLY,EDGE)
209 XV2 = X_VERTEX(POLY,EDGE + 1)
210
211 YV1 = Y_VERTEX(POLY,EDGE)
212 YV2 = Y_VERTEX(POLY,EDGE + 1)
213
214
215 D1 = DSQRT((x1 - XV1)**2 + (x2 - YV1)**2)
216 D2 = DSQRT((x1 - XV2)**2 + (x2 - YV2)**2)
217 IF(DMIN1(D1,D2)<TOL_POLY) THEN
218 F_POL = ZERO
219 BCID = BC_ID_P(POLY,EDGE)
220 RETURN
221 ENDIF
222
223 IF(DABS(YV2 - YV1) < TOL_POLY) THEN
224 test1 = (DABS(x2 - YV1) < TOL_POLY)
225 x_west = DMIN1(XV1,XV2)
226 x_east = DMAX1(XV1,XV2)
227 test2 = (x_west <= x1).AND.( x1 <= x_east)
228 IF (test1.and.test2) THEN
229 F_POL = ZERO
230 BCID = BC_ID_P(POLY,EDGE)
231 RETURN
232 ENDIF
233 ELSE
234 y_south = DMIN1(YV1,YV2)
235 y_north = DMAX1(YV1,YV2)
236 test1 = (y_south <= x2).AND.( x2 <= y_north)
237 IF(test1) THEN
238 slope = (XV2 - XV1) / (YV2 - YV1)
239 x_star = XV1 + slope * (x2 - YV1)
240
241 IF(DABS(slope)>max_slope) THEN
242 TOL_xstar = max_slope * TOL_POLY
243 ELSE
244 TOL_xstar = TOL_POLY
245 ENDIF
246
247 IF (x_star < x1 - TOL_xstar) THEN
248 COUNTER =COUNTER + 1
249 ELSEIF(DABS(x_star - x1) <= TOL_xstar) THEN
250 F_POL = ZERO
251 BCID = BC_ID_P(POLY,EDGE)
252 RETURN
253 ENDIF
254 ENDIF
255 ENDIF
256 ENDDO
257
258
259
260 IF((COUNTER/2) * 2 == COUNTER) THEN
261 F_POLY(POLY) = -POLY_SIGN(POLY)
262 ELSE
263 F_POLY(POLY) = POLY_SIGN(POLY)
264 ENDIF
265
266 ELSE
267 F_POLY(POLY) = -POLY_SIGN(POLY)
268 ENDIF
269
270 ENDDO
271
272
273 f_pol = MAXVAL(F_POLY(1:N_POLYGON))
274
275
276 CLIP_FLAG = .TRUE.
277
278
279 RETURN
280
281
282 END SUBROUTINE EVAL_POLY_FCT
283