psmile_overlap_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PSMILe_overlap_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_overlap_dble ( n_corners_src, n_corners_tgt, src, tgt, overlap )
00012 !
00013 ! !USES:
00014 !
00015         use PSMILe
00016 !
00017         implicit none
00018 !
00019 ! !INPUT PARAMETERS:
00020 !
00021         Integer, Intent(In)           :: n_corners_src
00022         Integer, Intent(In)           :: n_corners_tgt
00023         Type (line_dble), Intent(In)  :: src(n_corners_src)
00024         Type (line_dble), Intent(In)  :: tgt(n_corners_tgt)
00025 !
00026 ! !OUTPUT PARAMETERS:
00027 !
00028         Logical, Intent(Out)          :: overlap
00029 !
00030 ! !LOCAL VARIABLES
00031 !
00032         Integer                       :: i, j
00033 
00034         Logical                       :: intersect_dble
00035         Logical                       :: inside_dble
00036         Logical                       :: exact_dble
00037 !
00038 ! !DESCRIPTION:
00039 !
00040 ! http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/index.html
00041 ! http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
00042 ! http://en.wikipedia.org/wiki/Point_in_polygon
00043 !
00044 ! !REVISION HISTORY:
00045 !
00046 !   Date      Programmer   Description
00047 ! ----------  ----------   -----------
00048 ! 07.12.04    R. Redler    created
00049 !
00050 !EOP
00051 !----------------------------------------------------------------------
00052 !
00053 ! $Id: $
00054 ! $Author: $
00055 !
00056    Character(len=len_cvs_string), save :: mycvs = 
00057        '$Id: $'
00058 !
00059 !----------------------------------------------------------------------
00060 !
00061 !  Initialization
00062 !
00063         overlap = .false.
00064 !
00065 !  Test intersection of two edges
00066 !
00067         do j = 1, n_corners_tgt
00068            do i = 1, n_corners_src
00069               overlap = intersect_dble ( src(i), tgt(j) )
00070               if ( overlap ) return
00071            enddo
00072         enddo
00073 !
00074 !  Test whether one point is inside the cell
00075 !
00076         overlap = inside_dble ( n_corners_src, n_corners_tgt, src, tgt )
00077         if ( overlap ) return
00078 !
00079 !  Test whether source and target cell have an exact match
00080 !
00081         if ( n_corners_src == n_corners_tgt ) &
00082            overlap = exact_dble ( n_corners_src, src, tgt )
00083 
00084       end subroutine psmile_overlap_dble
00085 
00086 
00087 
00088       logical function intersect_dble ( line1, line2 )
00089 
00090          use PSMILe, only : shlon, line_dble, zero, one
00091          implicit none
00092 
00093          type(line_dble) :: line1
00094          type(line_dble) :: line2
00095 
00096          Double Precision, Parameter :: tol = 1.0d-12
00097 
00098          Double Precision :: determinant
00099          Double Precision :: numer1
00100          Double Precision :: numer2
00101 
00102          integer          :: i
00103 
00104          ! For calculating the intersection of two lines see e.g.
00105          ! http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/index.html
00106 
00107          intersect_dble = .false.
00108 
00109          determinant = ( line2%p2%y - line2%p1%y ) * ( line1%p2%x - line1%p1%x ) - &
00110                        ( line2%p2%x - line2%p1%x ) * ( line1%p2%y - line1%p1%y )
00111 
00112          if ( abs(determinant) < tol ) return ! lines are parallel
00113 
00114          ! We add or subtract 360.0 deg from the x-coordinates of
00115          ! the two points belonging to line 1 in addition to the
00116          ! standard case
00117 
00118          do i = 0, 2
00119 
00120             numer1 = ( line2%p2%x - line2%p1%x                ) * &
00121                      ( line1%p1%y - line2%p1%y                ) - &
00122                      ( line2%p2%y - line2%p1%y                ) * &
00123                      ( line1%p1%x - line2%p1%x +     shlon(i) )
00124 
00125             numer2 = ( line1%p2%x - line1%p1%x                ) * &
00126                      ( line1%p1%y - line2%p1%y                ) - &
00127                      ( line1%p2%y - line1%p1%y                ) * &
00128                      ( line1%p1%x - line2%p1%x +     shlon(i) )
00129 
00130             if ( ( numer1/determinant >= zero .and. numer1/determinant <= one ) .and. &
00131                  ( numer2/determinant >= zero .and. numer2/determinant <= one ) ) then
00132                    intersect_dble = .true.
00133 
00134                    return
00135             endif
00136 
00137          enddo
00138 
00139       end function intersect_dble
00140 
00141 
00142 
00143       logical function inside_dble ( nbr_corners_src, nbr_corners_tgt, src, tgt )
00144 
00145          use PSMILe, only : shlon, line_dble
00146          implicit none
00147 
00148          Integer, Intent(In)           :: nbr_corners_src
00149          Integer, Intent(In)           :: nbr_corners_tgt
00150          Type (line_dble), Intent(In)  :: src(nbr_corners_src)
00151          Type (line_dble), Intent(In)  :: tgt(nbr_corners_tgt)
00152 
00153          Double Precision :: nominator, denominator, xinters
00154 
00155          integer          :: i, k, n, counter
00156 
00157          ! Here we use a ray casting algorithm, see e.g.
00158          ! http://en.wikipedia.org/wiki/Point_in_polygon
00159          ! http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
00160 
00161          inside_dble = .false.
00162 
00163          ! We add or subtract 360.0 deg from the x-coordinates of the element
00164          ! for which we do the investigation in the n-loop
00165 
00166          do n = 0, 2
00167 
00168             ! First we check whether any of the target corners falls
00169             ! into the source element. In case of dealing with parallel
00170             ! edges between source and target cells this test is more
00171             ! robust than the above (testing for intersection between any
00172             ! two edges).
00173             
00174             do k = 1, nbr_corners_tgt
00175 
00176                counter = 0
00177 
00178                do i = 1, nbr_corners_src
00179 #if 0
00180                   if ( (src(i)%p1%y <= tgt(k)%p1%y  .and. &
00181                         tgt(k)%p1%y <  src(i)%p2%y) .or.  &
00182                        (src(i)%p2%y <= tgt(k)%p1%y  .and. &
00183                         tgt(k)%p1%y <  src(i)%p1%y) ) then
00184                      nominator   = (src(i)%p2%x-src(i)%p1%x) * (tgt(k)%p1%y-src(i)%p1%y)
00185                      denominator =  src(i)%p2%y-src(i)%p1%y
00186                      xinters     = nominator/denominator + src(i)%p1%x
00187                      if ( tgt(k)%p1%x < xinters ) counter = counter + 1
00188                   endif
00189 #else
00190                   if ( tgt(k)%p1%y > min(src(i)%p1%y,src(i)%p2%y) ) then
00191                      if ( tgt(k)%p1%y <= max(src(i)%p1%y,src(i)%p2%y) ) then
00192                         if ( tgt(k)%p1%x + shlon(n) <= max(src(i)%p1%x,src(i)%p2%x) ) then 
00193                            if ( src(i)%p1%y /= src(i)%p2%y ) then
00194                               nominator   = (tgt(k)%p1%y-src(i)%p1%y) * (src(i)%p2%x-src(i)%p1%x)
00195                               denominator = (src(i)%p2%y-src(i)%p1%y)
00196                               xinters     = nominator / denominator + src(i)%p1%x
00197                               if ( src(i)%p1%x == src(i)%p2%x .or. tgt(k)%p1%y <= xinters ) counter = counter + 1
00198                            endif
00199                         endif
00200                      endif
00201                   endif
00202 #endif
00203                enddo
00204 
00205                if ( mod(counter,2) /= 0 ) then
00206                   inside_dble = .true.
00207                   return
00208                endif
00209 
00210             enddo
00211 
00212             ! Second we check whether any of the source corners falls
00213             ! into the target element.
00214             
00215             do k = 1, nbr_corners_src
00216 
00217                counter = 0
00218 
00219                do i = 1, nbr_corners_tgt
00220 #if 0
00221                   if ( (tgt(i)%p1%y <= src(k)%p1%y  .and. &
00222                         src(k)%p1%y <  tgt(i)%p2%y) .or.  &
00223                        (tgt(i)%p2%y <= src(k)%p1%y  .and. &
00224                         src(k)%p1%y <  tgt(i)%p1%y) ) then
00225                      nominator   = (tgt(i)%p2%x-tgt(i)%p1%x) * (src(k)%p1%y-tgt(i)%p1%y)
00226                      denominator =  tgt(i)%p2%y-tgt(i)%p1%y
00227                      xinters     = nominator/denominator + tgt(i)%p1%x
00228                      if ( src(k)%p1%x < xinters ) counter = counter + 1
00229                   endif
00230 #else
00231                   if ( src(k)%p1%y > min(tgt(i)%p1%y,tgt(i)%p2%y) ) then
00232                      if ( src(k)%p1%y <= max(tgt(i)%p1%y,tgt(i)%p2%y) ) then
00233                         if ( src(k)%p1%x + shlon(n) <= max(tgt(i)%p1%x,tgt(i)%p2%x) ) then 
00234                            if ( tgt(i)%p1%y /= tgt(i)%p2%y ) then
00235                               nominator   = (src(k)%p1%y-tgt(i)%p1%y) * (tgt(i)%p2%x-tgt(i)%p1%x)
00236                               denominator = (tgt(i)%p2%y-tgt(i)%p1%y) 
00237                               xinters     = nominator / denominator +  tgt(i)%p1%x
00238                               if ( tgt(i)%p1%x == tgt(i)%p2%x .or. src(k)%p1%x <= xinters ) counter = counter + 1
00239                            endif
00240                         endif
00241                      endif
00242                   endif
00243 #endif
00244                enddo
00245 
00246                if ( mod(counter,2) /= 0 ) then
00247                   inside_dble = .true.
00248                   return
00249                endif
00250 
00251             enddo
00252 
00253          enddo
00254 
00255       end function inside_dble
00256 
00257 
00258 
00259       logical function exact_dble ( nbr_corners, src, tgt )
00260 
00261          use PSMILe, only : shlon, line_dble
00262          implicit none
00263 
00264          Integer, Intent(In)           :: nbr_corners
00265          Type (line_dble), Intent(In)  :: src(nbr_corners)
00266          Type (line_dble), Intent(In)  :: tgt(nbr_corners)
00267 
00268          Double Precision :: tol = 1.0d-12
00269 
00270          integer          :: i, j, k, n
00271 
00272          exact_dble = .false.
00273 
00274          ! We add or subtract 360.0 deg from the x-coordinates of
00275          ! the two points belonging the source point 1 in addition
00276          ! to the standard case.
00277 
00278          do k = 0, 2
00279             n = 0
00280             do i = 1, nbr_corners
00281                do j = 1, nbr_corners
00282                   if ( abs(src(i)%p1%x + shlon(k) - tgt(j)%p1%x) < tol .and. &
00283                        abs(src(i)%p1%y            - tgt(j)%p1%y) < tol ) then
00284                      n = n + 1
00285                      exit
00286                   endif
00287                enddo
00288             enddo
00289             exact_dble = ( n > 2 ) ! this may not detect all cases
00290             if ( exact_dble ) return
00291          enddo
00292 
00293       end function exact_dble