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 + shlon(n) < 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. &
00198                                    tgt(k)%p1%x + shlon(n) <= xinters ) counter = counter + 1
00199                            endif
00200                         endif
00201                      endif
00202                   endif
00203 #endif
00204                enddo
00205 
00206                if ( mod(counter,2) /= 0 ) then
00207                   inside_dble = .true.
00208                   return
00209                endif
00210 
00211             enddo
00212 
00213             ! Second we check whether any of the source corners falls
00214             ! into the target element.
00215             
00216             do k = 1, nbr_corners_src
00217 
00218                counter = 0
00219 
00220                do i = 1, nbr_corners_tgt
00221 #if 0
00222                   if ( (tgt(i)%p1%y <= src(k)%p1%y  .and. &
00223                         src(k)%p1%y <  tgt(i)%p2%y) .or.  &
00224                        (tgt(i)%p2%y <= src(k)%p1%y  .and. &
00225                         src(k)%p1%y <  tgt(i)%p1%y) ) then
00226                      nominator   = (tgt(i)%p2%x-tgt(i)%p1%x) * (src(k)%p1%y-tgt(i)%p1%y)
00227                      denominator =  tgt(i)%p2%y-tgt(i)%p1%y
00228                      xinters     = nominator/denominator + tgt(i)%p1%x
00229                      if ( src(k)%p1%x + shlon(n) < xinters ) counter = counter + 1
00230                   endif
00231 #else
00232                   if ( src(k)%p1%y > min(tgt(i)%p1%y,tgt(i)%p2%y) ) then
00233                      if ( src(k)%p1%y <= max(tgt(i)%p1%y,tgt(i)%p2%y) ) then
00234                         if ( src(k)%p1%x + shlon(n) <= max(tgt(i)%p1%x,tgt(i)%p2%x) ) then 
00235                            if ( tgt(i)%p1%y /= tgt(i)%p2%y ) then
00236                               nominator   = (src(k)%p1%y-tgt(i)%p1%y) * (tgt(i)%p2%x-tgt(i)%p1%x)
00237                               denominator = (tgt(i)%p2%y-tgt(i)%p1%y)
00238                               xinters     = nominator / denominator + tgt(i)%p1%x
00239                               if ( tgt(i)%p1%x == tgt(i)%p2%x .or. &
00240                                    src(k)%p1%x + shlon(n) <= xinters ) counter = counter + 1
00241                            endif
00242                         endif
00243                      endif
00244                   endif
00245 #endif
00246                enddo
00247 
00248                if ( mod(counter,2) /= 0 ) then
00249                   inside_dble = .true.
00250                   return
00251                endif
00252 
00253             enddo
00254 
00255          enddo
00256 
00257       end function inside_dble
00258 
00259 
00260 
00261       logical function exact_dble ( nbr_corners, src, tgt )
00262 
00263          use PSMILe, only : shlon, line_dble
00264          implicit none
00265 
00266          Integer, Intent(In)           :: nbr_corners
00267          Type (line_dble), Intent(In)  :: src(nbr_corners)
00268          Type (line_dble), Intent(In)  :: tgt(nbr_corners)
00269 
00270          Double Precision :: tol = 1.0d-12
00271 
00272          integer          :: i, j, k, n
00273 
00274          exact_dble = .false.
00275 
00276          ! We add or subtract 360.0 deg from the x-coordinates of
00277          ! the two points belonging the source point 1 in addition
00278          ! to the standard case.
00279 
00280          do k = 0, 2
00281             n = 0
00282             do i = 1, nbr_corners
00283                do j = 1, nbr_corners
00284                   if ( abs(src(i)%p1%x + shlon(k) - tgt(j)%p1%x) < tol .and. &
00285                        abs(src(i)%p1%y            - tgt(j)%p1%y) < tol ) then
00286                      n = n + 1
00287                      exit
00288                   endif
00289                enddo
00290             enddo
00291             exact_dble = ( n > 2 ) ! this may not detect all cases
00292             if ( exact_dble ) return
00293          enddo
00294 
00295       end function exact_dble