get_poly_data.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_POLY_DATA                                          C
!  Purpose: reads polygon(s) coordinates from poly.dat (2D only)       C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 30-JAN-09  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!  Revision Number:                                                    C
!  Purpose:                                                            C
!  Author:                                            Date: dd-mmm-yy  C
!  Reviewer:                                          Date: dd-mmm-yy  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced:                                               C
!                                                                      C
!  Variables modified:                                                 C
!                                                                      C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
!
      SUBROUTINE get_poly_data

!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------

      USE compar
      USE exit, only: mfix_exit
      USE fldvar
      USE funits
      USE mpi_utility
      USE param
      USE param1
      USE physprop
      USE polygon
      USE progress_bar
      USE run
      USE rxns
      USE scalars
      IMPLICIT NONE

      INTEGER :: POLY,V,NN,NSKIP
      LOGICAL :: PRESENT

      WRITE(*,2000) 'READING polygon geometry from poly.dat...'

      INQUIRE(file='poly.dat',exist=present)
      IF(.NOT.present) THEN
         IF(mype == pe_io) THEN
            WRITE(*,"('(PE ',I3,'): input data file, ',A11,' is missing: run aborted')") &
            mype,'poly.dat'
         ENDIF
         CALL mfix_exit(mype)
      ENDIF
!
!
!     OPEN poly.dat ASCII FILE
!
      OPEN(convert='BIG_ENDIAN',unit=333, file='poly.dat', status='OLD', err=910)

      nskip = 13
      DO nn=1,nskip
         READ(333,*,err=920,end=930)   ! Skip first NSKIP comment lines
      ENDDO

      READ(333,*,err=920,end=930) n_polygon

      DO poly = 1, n_polygon
         READ(333,*,err=920,end=930) n_vertex(poly), poly_sign(poly)
         IF(dabs(poly_sign(poly))/=one) THEN
            WRITE(*,*)'ERROR WHILE READIND poly.dat:'
            WRITE(*,*)'POLYGON SIGN MUST BE +1.0 or -1.0.'
            CALL mfix_exit(mype)
         ENDIF
         DO v = 1,n_vertex(poly)
            READ(333,*,err=920,end=930) x_vertex(poly,v),y_vertex(poly,v),bc_id_p(poly,v)
         ENDDO
         x_vertex(poly,n_vertex(poly) + 1) = x_vertex(poly, 1)  ! Copy first vertex as last vertex to close the polygon
         y_vertex(poly,n_vertex(poly) + 1) = y_vertex(poly, 1)
      ENDDO

      CLOSE(333)

      WRITE(*,2010) 'Polygon geometry successfully read for ', n_polygon, ' polygon(s).'

      RETURN

!======================================================================
!     HERE IF AN ERROR OCCURED OPENNING/READING THE FILE
!======================================================================
!
 910  CONTINUE
      WRITE (*, 1500)
      CALL mfix_exit(mype)
 920  CONTINUE
      WRITE (*, 1600)
      CALL mfix_exit(mype)
 930  CONTINUE
      WRITE (*, 1700)
      CALL mfix_exit(mype)
!
 1500 FORMAT(/1x,70('*')//' From: GET_POLY_DATA',/' Message: ',&
      'Unable to open poly.dat file',/1x,70('*')/)
 1600 FORMAT(/1x,70('*')//' From: GET_POLY_DATA',/' Message: ',&
      'Error while reading poly.dat file',/1x,70('*')/)
 1700 FORMAT(/1x,70('*')//' From: GET_POLY_DATA',/' Message: ',&
      'End of file reached while reading poly.dat file',/1x,70('*')/)
 2000 FORMAT(1x,a)
 2010 FORMAT(1x,a,i4,a)

      END SUBROUTINE get_poly_data


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: EVAL_POLY_FCT                                          C
!  Purpose: Evaluates a user-defined function                          C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 21-Feb-08  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number #                                  Date: ##-###-##  C
!  Author: #                                                           C
!  Purpose: #                                                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE eval_poly_fct(x1,x2,x3,Q,f_pol,CLIP_FLAG,BCID)

      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE fldvar
      USE quadric
      USE cutcell
      USE polygon

      IMPLICIT NONE
      DOUBLE PRECISION x1,x2,x3
      DOUBLE PRECISION f_pol
      DOUBLE PRECISION, DIMENSION(DIM_POLYGON) :: F_POLY
      INTEGER :: Q
      LOGICAL :: CLIP_X,CLIP_Y,CLIP_FLAG
      LOGICAL :: test1,test2

      INTEGER :: POLY,N_EDGE,COUNTER,EDGE
      DOUBLE PRECISION :: X_VERTEX_MIN,X_VERTEX_MAX,Y_VERTEX_MIN,Y_VERTEX_MAX
      DOUBLE PRECISION :: XV1,XV2,YV1,YV2
      DOUBLE PRECISION :: x_west,x_east,y_north,y_south,slope,x_star,D1,D2
      DOUBLE PRECISION :: max_slope = 200.0d0
      DOUBLE PRECISION :: TOL_xstar
      INTEGER :: BCID

!======================================================================
! This subroutine checks whether a point P(x1,x2,x3) lies inside any
! of the N_POLGON polygon(s), using the Shimrat's algorithm:
! 1) Draw a horizontal line from point P to -infinity (here taken as
!    the minimum x-location of the polygon vertices)
! 2) Count the number of intersection(s) with all edges (COUNTER)
! 3) If COUNTER is even, the point is outside
!    If COUNTER is odd,  the point is inside
! Once the point's location is determined, a value is assigned to f_pol:
!    fpol = 0 if point lies on top of one of the polygon's edges
!                (within tolerance TOL_POLY)
!    fpol = POLY_SIGN(POLY) if point lies inside the polygon
!    fpol = - POLY_SIGN(POLY) if point lies outside the polygon
!    where
!           POLY_SIGN(POLY) = -1.0 means the interior of the polygon
!                                  is part of the computational domain
!           POLY_SIGN(POLY) = +1.0 means the interior of the polygon
!                                  is excluded from the computational domain
!======================================================================

      IF(n_polygon < 1) RETURN

      bcid = -1

      DO poly = 1,n_polygon

         n_edge   = n_vertex(poly)

         x_vertex_min = minval(x_vertex(poly,:)) - 2.0d0 * tol_poly
         x_vertex_max = maxval(x_vertex(poly,:)) + 2.0d0 * tol_poly

         y_vertex_min = minval(y_vertex(poly,:)) - 2.0d0 * tol_poly
         y_vertex_max = maxval(y_vertex(poly,:)) + 2.0d0 * tol_poly

         clip_x = (x_vertex_min <= x1).AND.( x1 <= x_vertex_max)
         clip_y = (y_vertex_min <= x2).AND.( x2 <= y_vertex_max)


         IF(clip_x.AND.clip_y) THEN

            counter = 0

            DO edge = 1, n_edge

               xv1 = x_vertex(poly,edge)
               xv2 = x_vertex(poly,edge + 1)

               yv1 = y_vertex(poly,edge)
               yv2 = y_vertex(poly,edge + 1)


               d1 = dsqrt((x1 - xv1)**2 + (x2 - yv1)**2)
               d2 = dsqrt((x1 - xv2)**2 + (x2 - yv2)**2)
               IF(dmin1(d1,d2)<tol_poly) THEN
                  f_pol = zero
                  bcid = bc_id_p(poly,edge)
                  RETURN
               ENDIF

               IF(dabs(yv2 - yv1) < tol_poly) THEN
                  test1 = (dabs(x2 - yv1) < tol_poly)
                  x_west = dmin1(xv1,xv2)
                  x_east = dmax1(xv1,xv2)
                  test2 = (x_west <= x1).AND.( x1 <= x_east)
                  IF (test1.and.test2) THEN
                     f_pol = zero
                     bcid = bc_id_p(poly,edge)
                     RETURN
                  ENDIF
               ELSE
                  y_south = dmin1(yv1,yv2)
                  y_north = dmax1(yv1,yv2)
                  test1 = (y_south <= x2).AND.( x2 <= y_north)
                  IF(test1) THEN
                     slope = (xv2 - xv1) / (yv2 - yv1)
                     x_star = xv1 + slope * (x2 - yv1)

                     IF(dabs(slope)>max_slope) THEN
                        tol_xstar = max_slope * tol_poly
                     ELSE
                        tol_xstar = tol_poly
                     ENDIF

                     IF (x_star < x1 - tol_xstar) THEN
                        counter =counter + 1
                     ELSEIF(dabs(x_star - x1) <= tol_xstar) THEN
                        f_pol = zero
                        bcid = bc_id_p(poly,edge)
                        RETURN
                     ENDIF
                  ENDIF
               ENDIF
            ENDDO

            IF((counter/2) * 2 == counter) THEN
               f_poly(poly) = -poly_sign(poly)    ! even ==> outside
            ELSE
               f_poly(poly) = poly_sign(poly)     ! odd  ==> inside
            ENDIF

         ELSE
            f_poly(poly) = -poly_sign(poly)       ! outside
         ENDIF

      ENDDO

      f_pol = maxval(f_poly(1:n_polygon))

      clip_flag = .true.

      RETURN

      END SUBROUTINE eval_poly_fct