File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/get_poly_data.f

1     
2     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
3     !                                                                      C
4     !  Module name: GET_POLY_DATA                                          C
5     !  Purpose: reads polygon(s) coordinates from poly.dat (2D only)       C
6     !                                                                      C
7     !  Author: Jeff Dietiker                              Date: 30-JAN-09  C
8     !  Reviewer:                                          Date: **-***-**  C
9     !                                                                      C
10     !  Revision Number:                                                    C
11     !  Purpose:                                                            C
12     !  Author:                                            Date: dd-mmm-yy  C
13     !  Reviewer:                                          Date: dd-mmm-yy  C
14     !                                                                      C
15     !  Literature/Document References:                                     C
16     !                                                                      C
17     !  Variables referenced:                                               C
18     !                                                                      C
19     !  Variables modified:                                                 C
20     !                                                                      C
21     !  Local variables:                                                    C
22     !                                                                      C
23     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
24     !
25           SUBROUTINE GET_POLY_DATA
26     
27     !-----------------------------------------------
28     !   M o d u l e s
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     !     OPEN poly.dat ASCII FILE
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)   ! Skip first NSKIP comment lines
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)  ! Copy first vertex as last vertex to close the polygon
83              Y_VERTEX(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     !     HERE IF AN ERROR OCCURED OPENNING/READING THE FILE
94     !======================================================================
95     !
96      910  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      1500 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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
119     !                                                                      C
120     !  Module name: EVAL_POLY_FCT                                          C
121     !  Purpose: Evaluates a user-defined function                          C
122     !                                                                      C
123     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
124     !  Reviewer:                                          Date:            C
125     !                                                                      C
126     !  Revision Number #                                  Date: ##-###-##  C
127     !  Author: #                                                           C
128     !  Purpose: #                                                          C
129     !                                                                      C
130     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
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     ! This subroutine checks whether a point P(x1,x2,x3) lies inside any
166     ! of the N_POLGON polygon(s), using the Shimrat's algorithm:
167     ! 1) Draw a horizontal line from point P to -infinity (here taken as
168     !    the minimum x-location of the polygon vertices)
169     ! 2) Count the number of intersection(s) with all edges (COUNTER)
170     ! 3) If COUNTER is even, the point is outside
171     !    If COUNTER is odd,  the point is inside
172     ! Once the point's location is determined, a value is assigned to f_pol:
173     !    fpol = 0 if point lies on top of one of the polygon's edges
174     !                (within tolerance TOL_POLY)
175     !    fpol = POLY_SIGN(POLY) if point lies inside the polygon
176     !    fpol = - POLY_SIGN(POLY) if point lies outside the polygon
177     !    where
178     !           POLY_SIGN(POLY) = -1.0 means the interior of the polygon
179     !                                  is part of the computational domain
180     !           POLY_SIGN(POLY) = +1.0 means the interior of the polygon
181     !                                  is excluded from the computational domain
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)    ! even ==> outside
262                 ELSE
263                    F_POLY(POLY) = POLY_SIGN(POLY)     ! odd  ==> inside
264                 ENDIF
265     
266              ELSE
267                 F_POLY(POLY) = -POLY_SIGN(POLY)       ! outside
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