DOUG 0.2

Polygon.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 module Polygon_class
00023 
00024   use DOUG_utils
00025   use Points2D_class
00026   use globals
00027 
00028   implicit none
00029 
00030 #include<doug_config.h>
00031 
00032 !!$  include 'DOUG_utils.f90'
00033 !!$  include 'globals.f90'
00034 !!$  include 'Points2D.f90'
00035 
00036   type Polygon
00037      type(Points2D) :: pts      ! Points in 2D forming a polygon
00038 
00039      logical :: closed=.false.
00040      logical :: sorted=.false.
00041   end type Polygon
00042 
00043   public :: &
00044        Polygon_New,       &
00045        Polygon_Destroy,   &
00046        Polygon_Init,      &
00047        Polygon_getSize,   &
00048        Polygon_getX,      &
00049        Polygon_getY,      &
00050        Polygon_sortVerts, &
00051        Polygon_isSorted,  &
00052        Polygon_Close,     &
00053        Polygon_isClosed,  &
00054        Polygon_Area,      &
00055        Polygon_Centroid,  &
00056        Polygon_pl2D_Plot, &
00057        Polygon_Print
00058 
00059   !private ::
00060 
00061 contains
00062 
00063 
00064   !--------------------------------
00065   ! Class constructor
00066   !--------------------------------
00067   function Polygon_New(n) result(p)
00068 
00069     integer        :: n
00070     type(Polygon)  :: p
00071 
00072 
00073     p%pts = Points2D_New(n)
00074 
00075     p%pts%x = 0
00076     p%pts%y = 0
00077 
00078   end function Polygon_New
00079 
00080 
00081   !----------------------------
00082   ! Class destructor
00083   !----------------------------
00084   subroutine Polygon_Destroy(p)
00085 
00086     type(Polygon) :: p
00087 
00088 
00089     call Points2D_Destroy(p%pts)
00090 
00091     p%closed = .false.
00092     p%sorted = .false.
00093 
00094   end subroutine Polygon_Destroy
00095 
00096 
00097   !-------------------------------------------------
00098   ! Initialise polygon
00099   !-------------------------------------------------
00100   subroutine Polygon_Init(p, xc, yc, sorted, closed)
00101 
00102     type(Polygon),               intent(in out) :: p
00103     real(kind=xyzk), dimension(:), intent(in)   :: xc, yc
00104     logical, intent(in), optional               :: sorted
00105     logical, intent(in), optional               :: closed
00106 
00107 
00108     call Polygon_Destroy(p)
00109 
00110     p%pts = Points2D_newFill(xc, yc)
00111 
00112     if (present(sorted)) p%sorted = sorted
00113     if (present(closed)) p%closed = closed
00114 
00115   end subroutine Polygon_Init
00116 
00117 
00118   !---------------------------------------
00119   ! Get number of points forming a polygon
00120   !---------------------------------------
00121   function Polygon_getSize(p) result(size)
00122 
00123     type(Polygon), intent(in) :: p
00124     integer                   :: size
00125 
00126     size = p%pts%np
00127 
00128   end function Polygon_getSize
00129 
00130 
00131   !---------------------------------
00132   ! Return x-coordinates
00133   !---------------------------------
00134   function Polygon_getX(p) result(x)
00135     implicit none
00136 
00137     type(Polygon), intent(in)              :: p
00138     real(kind=xyzk), dimension(:), pointer :: x
00139 
00140     allocate(x(Polygon_getSize(p)))
00141 
00142     x = p%pts%x
00143 
00144   end function Polygon_getX
00145 
00146 
00147   !---------------------------------
00148   ! Return y-coordinates
00149   !---------------------------------
00150   function Polygon_getY(p) result(y)
00151     implicit none
00152 
00153     type(Polygon), intent(in)              :: p
00154     real(kind=xyzk), dimension(:), pointer :: y
00155 
00156     allocate(y(Polygon_getSize(p)))
00157 
00158     y = p%pts%y
00159 
00160   end function Polygon_getY
00161 
00162 
00163   !--------------------------------------------------------
00164   ! Sort vertices in polygon (via reodering by coordinates)
00165   ! to make it convex (not self intersecting). Applicable
00166   ! in case of points forming an element.
00167   ! NB:
00168   ! Actually, elements are allready "convex", but their
00169   ! nodes must not necessarily be represented in a
00170   ! sorted way to form convex hulls.
00171   !--------------------------------------------------------
00172   subroutine Polygon_sortVerts(p)
00173 
00174     use globals, only : stream, D_PI25DT, D_PI2
00175 
00176     implicit none
00177 
00178     type(Polygon), intent(in out)  :: p ! Polygon representing an element
00179 
00180     type(Points2D)                 :: p1
00181     integer                        :: p1_ind
00182     real(kind=xyzk)                :: x_min, y_min, temp, tmpx, tmpy
00183     real(kind=xyzk), dimension(:), allocatable  :: xc, yc
00184     real(kind=xyzk), dimension(:), allocatable  :: angles
00185     integer                        :: count, i, j
00186     integer                        :: n
00187     integer, dimension(:), pointer :: min_inds
00188     character*10                   :: buf10
00189 
00190     real(kind=xyzk)                :: Degr = 180.0_xyzk / D_PI25DT
00191 
00192 
00193     if (Polygon_isClosed(p)) then
00194        call DOUG_abort('[Polygon_sortVerts] : can not sort closed polygons yet - TODO.')
00195     end if
00196 
00197     if (Polygon_isSorted(p)) return
00198 
00199     n = Polygon_getSize(p)
00200 
00201     allocate(xc(n), yc(n))
00202     xc = Polygon_getX(p)
00203     yc = Polygon_getY(p)
00204 
00205     ! Position polygon in 1st quadrant :
00206     ! (with min(y)=0.0, so laying on X-axis and min(x)=0.0 - laying on Y-axis)
00207     yc = yc - minval(yc)
00208     xc = xc - minval(xc)
00209 
00210 
00211     ! Find extremal point to work with.
00212     !
00213     ! The one with smallest y coordinate. If there are many -
00214     ! choose the one with smallest x coordinate among of them.
00215     !
00216 
00217     ! Get y minimum.
00218     ! NB: Due to roundoff errors on some machines
00219     ! it must not necessarily be equal to 0.0_xyzk
00220     y_min = minval(yc) ! y_min = 0.0_xyzk
00221     x_min = minval(xc) ! x_min = 0.0_xyzk
00222 
00223     ! Count number of min. points. There could be more than one.
00224     allocate(min_inds(n)); min_inds = 0
00225     count = 0
00226     do i = 1,n
00227        if (y_min == yc(i)) then
00228           count = count + 1
00229           min_inds(count) = i
00230        end if
00231     end do
00232  !   write(stream, *) 'Number of y_min points is ', count
00233 
00234     if (count == 1) then
00235        ! Only one minimum
00236        p1_ind = min_inds(count)
00237        p1 = Points2D_newFill(xc(p1_ind:p1_ind), yc(p1_ind:p1_ind))
00238     else if (count > 1) then
00239        ! We are having more than one minimum
00240  !      write(stream, *) 'We are having multple ''Y'' minima in : ', min_inds(:count)
00241  !      write(stream, *) 'with ''X'' values : ', xc(min_inds(:count))
00242        ! Find the one with the smalest x coordinate
00243        x_min = minval(xc(min_inds(:count)))
00244  !      write(stream, *) 'smallest is ', x_min
00245        ! Get its index
00246        do i = 1,count
00247           if (x_min == xc(min_inds(i))) then
00248              p1_ind = min_inds(i)
00249              p1 = Points2D_newFill(xc(p1_ind:p1_ind), yc(p1_ind:p1_ind))
00250              exit
00251           end if
00252        end do
00253 !       write(stream, *) 'its index is ', p1_ind
00254     else
00255        write(buf10, '(i10)') count
00256        call DOUG_abort('[Polygon_sortVerts] : some weird value of count = '//buf10//'.')
00257     end if
00258 
00259 !    write(stream, *) 'will work with the vertex ', p1%x, p1%y
00260 
00261 
00262     ! Find angles between X-axis and [p_1,p_i] lines, where i=2..n.
00263     allocate(angles(n))
00264     do i = 1,n
00265        angles(i) = Degr * atan((yc(i) - p1%y(1)) / (xc(i) - p1%x(1)))
00266 
00267 !       write(stream, *) 'angles(',i,') = ', angles(i)
00268 
00269        if (angles(i) < 0) then ! Obtuse angle
00270           angles(i) = 180.0_xyzk + angles(i)
00271 
00272 !          write(stream, *) 'angles_i<0 : after: angles(',i,') = ', angles(i)
00273 
00274        else if (isinf(angles(i))==1) then ! Infinity
00275 
00276           ! As tan(90) is equal to infinity, substitute it with
00277           ! the corresponding value in degrees.
00278           angles(i) = 90.0_xyzk
00279 
00280 !          write(stream, *) 'angles_i==Inf : after: angles(',i,') = ', angles(i)
00281 
00282        else if (isnan(angles(i))) then
00283 
00284           ! WARNING!
00285           !
00286           ! HACK to simplify sorting algorythm.
00287           ! Set p_1 point to be equal to -1.0. Ascending sorting will allways
00288           ! put it on the first place because all angle values are > zero.
00289 
00290           angles(i) = -1.0_xyzk
00291 
00292 !          write(stream, *) 'angles_i<FIRST : after: angles(',i,') = ', angles(i)
00293 
00294        end if
00295     end do
00296 
00297 !    write(stream, *) 'angles:', angles
00298 
00299 
00300     ! Sort points radially, using 'angles' and p_1 as the origin.
00301     !
00302     ! Algorithm - Insertion Sort:
00303 
00304     ! for i= n-1 down to 1 begin
00305     !   temp = x_i
00306     !   j = i+1
00307     !   while(j <= n and x_j < temp) begin
00308     !     x_(j-1) = x_j
00309     !     j=j+1
00310     !   end
00311     !   x_(j-1) = temp
00312     ! end
00313 
00314     do i = n-1,1,-1
00315        !write(stream, *) ' i = ', i
00316        temp = angles(i)
00317        tmpx = p%pts%x(i)
00318        tmpy = p%pts%y(i)
00319        j = i + 1
00320        do while (angles(j) < temp)
00321           angles(j-1) = angles(j)
00322           p%pts%x(j-1) = p%pts%x(j)
00323           p%pts%y(j-1) = p%pts%y(j)
00324           j = j + 1
00325           if (j > n) exit
00326        end do
00327        angles(j-1) = temp
00328        p%pts%x(j-1) = tmpx
00329        p%pts%y(j-1) = tmpy
00330        !write(stream, *) 'angles: ', angles
00331     end do
00332 !    write(stream, *) 'angles after sort: ', angles
00333 
00334 !    write(stream, *) 'p%pts%x:', p%pts%x
00335 !    write(stream, *) 'p%pts%y:', p%pts%y
00336 
00337     p%sorted = .true.
00338 
00339     deallocate(min_inds, xc, yc, angles)
00340     call Points2D_Destroy(p1)
00341 
00342   end subroutine Polygon_sortVerts
00343 
00344 
00345   !-----------------------------------------
00346   ! Is polygon sorted or not - test function
00347   !-----------------------------------------
00348   function Polygon_isSorted(p) result(res)
00349 
00350     type(Polygon), intent(in out) :: p
00351 
00352     logical :: res
00353 
00354     if (p%sorted.eqv.(.true.)) then
00355        res = .true.
00356     else
00357        res = .false.
00358     end if
00359   end function Polygon_isSorted
00360 
00361 
00362   !---------------------------------------
00363   ! Close set of points forming an element
00364   ! to make a closed polygon of it.
00365   !---------------------------------------
00366   subroutine Polygon_Close(p)
00367 
00368     implicit none
00369     type(Polygon), intent(in out) :: p
00370 
00371 
00372     if (Polygon_isClosed(p)) return
00373 
00374     call Points2D_resizePreserve(p%pts, p%pts%np+1)
00375 
00376     p%pts%x(p%pts%np) = p%pts%x(1)
00377     p%pts%y(p%pts%np) = p%pts%y(1)
00378 
00379     p%closed = .true.
00380 
00381   end subroutine Polygon_Close
00382 
00383 
00384   !-----------------------------------------
00385   ! Is polygon closed or not - test function
00386   !-----------------------------------------
00387   function Polygon_isClosed(p) result(res)
00388     implicit none
00389     type(Polygon), intent(in out) :: p
00390 
00391     logical :: res
00392 
00393     if (p%closed.eqv.(.true.)) then
00394        res = .true.
00395     else
00396        res = .false.
00397     end if
00398   end function Polygon_isClosed
00399 
00400 
00401   !-------------------------------------------
00402   ! Find area of a polygon
00403   !
00404   ! A = 1/2 * SUM[(Xi * Yi+1 - Xi+1 * Yi)]
00405   !
00406   ! NB:
00407   !   Polygon :
00408   !     * with no holes
00409   !     * must not be self intersecting
00410   !       (must be sorted to form convex hull)
00411   !-------------------------------------------
00412   function Polygon_Area(p) result(area)
00413 
00414     type(Polygon)   :: p
00415     real(kind=xyzk) :: area
00416 
00417     integer       :: n, i, j
00418 
00419     if (.not.Polygon_isSorted(p)) call Polygon_sortVerts(p)
00420 
00421     n = Polygon_getSize(p)
00422     if (Polygon_isClosed(p)) n = n - 1
00423 
00424     area = 0.0_xyzk
00425     do i = 1,n
00426        if (i == n) then
00427           j = 1
00428        else
00429           j = i+1
00430        end if
00431        area = area + (p%pts%x(i) * p%pts%y(j) - p%pts%y(i) * p%pts%x(j));
00432     end do
00433     area =  area / 2.0_xyzk;
00434 
00435     if (area < 0.0_xyzk) area = (-1.0_xyzk) * area
00436 
00437   end function Polygon_Area
00438 
00439 
00440   !--------------------------------------------------------
00441   ! Find center of mass of a polygon
00442   !
00443   !  X = SUM[(Xi + Xi+1) * (Xi * Yi+1 - Xi+1 * Yi)] / 6 / A
00444   !  Y = SUM[(Yi + Yi+1) * (Xi * Yi+1 - Xi+1 * Yi)] / 6 / A
00445   !--------------------------------------------------------
00446   function Polygon_Centroid(p) result(pcent)
00447 
00448     use globals, only : stream
00449 
00450     type(Polygon)  :: p
00451     type(Points2D) :: pcent
00452 
00453     real(kind=xyzk)  :: second_factor, Area
00454     integer          :: n, i, j
00455 
00456     if (.not.Polygon_isSorted(p)) call Polygon_sortVerts(p)
00457 
00458     n = Polygon_getSize(p)
00459     if (Polygon_isClosed(p)) n = n - 1
00460 
00461     pcent = Points2D_New(1)
00462     call Points2D_Init(pcent, (/0.0_xyzk/), (/0.0_xyzk/))
00463 
00464     do i = 1,n
00465        if (i == n) then
00466           j = 1
00467        else
00468           j = i+1
00469        end if
00470        second_factor = (p%pts%x(i) * p%pts%y(j) - p%pts%y(i) * p%pts%x(j))
00471        pcent%x = pcent%x + (p%pts%x(i) + p%pts%x(j)) * second_factor
00472        pcent%y = pcent%y + (p%pts%y(i) + p%pts%y(j)) * second_factor
00473     end do
00474 
00475     Area = Polygon_Area(p)
00476 
00477     ! Divide by 6 times the polygon's area
00478     pcent%x = pcent%x / 6.0_xyzk / Area
00479     pcent%y = pcent%y / 6.0_xyzk / Area
00480 
00481   end function Polygon_Centroid
00482   !=============================================
00483   !
00484   ! Plotting subroutines
00485   !
00486   !---------------------------------------------
00487   ! Plot polygon in 2D
00488   !---------------------------------------------
00489   subroutine Polygon_pl2D_Plot(p, INIT_CONT_END)
00490 
00491     use globals, only : stream, D_MSGLVL
00492 
00493     implicit none
00494 
00495     type(Polygon)                        :: p
00496     integer,      intent(in), optional   :: INIT_CONT_END
00497 
00498     real(kind=xyzk), dimension(:), pointer :: xc, yc
00499     real(kind=xyzk)                        :: xmin, xmax, ymin, ymax
00500     integer                              :: n
00501     character*3                          :: buf3
00502 
00503 #ifdef D_WANT_PLPLOT_YES
00504 
00505     if (D_MSGLVL > 5) then
00506        write(stream, *)
00507        write(stream, *) '[Polygon_pl2D_Plot] : Plotting 2D plygon.'
00508     end if
00509 
00510     if (.not.p%closed) then
00511        if (D_MSGLVL > 5) &
00512             write(stream, FMT='(a)', advance='no') '   Only closed polygons can be plotted. Closing plygon ... '
00513        call Polygon_Close(p)
00514        if (D_MSGLVL > 5) &
00515             write(stream, *) 'done'
00516     end if
00517 
00518     n = Polygon_getSize(p) ! Length of coordinate arrays
00519 
00520     xc => Polygon_getX(p)
00521     yc => Polygon_getY(p)
00522 
00523     if (.not.present(INIT_CONT_END).or.&
00524          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_INIT))) then
00525        xmin = minval(xc)
00526        xmax = maxval(xc)
00527        ymin = minval(yc)
00528        ymax = maxval(yc)
00529 
00530        xmin = xmin - abs(xmax)/10.0
00531        xmax = xmax + abs(xmax)/10.0
00532        ymin = ymin - abs(ymax)/10.0
00533        ymax = ymax + abs(ymax)/10.0
00534 
00535        call plsdev("tk")
00536        call plinit()
00537 
00538        call plenv (xmin, xmax, ymin, ymax, 0, 0);
00539 
00540        call plcol0(1) ! red
00541        write(buf3, '(i3)') n-1 ! Because we plot closed polygons
00542        call pllab( '(x)', '(y)', 'Polygon : '//buf3//' vertices' )
00543     end if
00544 
00545     ! Plot polygon
00546     call plcol0(3) ! green
00547     call plline(n, xc, yc)
00548 
00549     ! Plot vertices
00550     call plcol0(1) ! red
00551     call plssym(0.0d0, 2.0d0)
00552     call plpoin(n, xc, yc, 1)
00553 
00554     if (.not.present(INIT_CONT_END).or.&
00555          (present(INIT_CONT_END).and.(INIT_CONT_END == D_PLPLOT_END))) then
00556        call plend()
00557     end if
00558 
00559     nullify(xc, yc)
00560 
00561 #else
00562     write(stream, '(/a)') ' [Polygon_pl2D_Plot] : Compiled w/o plotting support!'
00563 #endif
00564 
00565 
00566   end subroutine Polygon_pl2D_Plot
00567   !=============================
00568   !
00569   ! I/O subroutines
00570   !
00571   !-----------------------------
00572   ! Print polygon for human ayes
00573   !-----------------------------
00574   subroutine Polygon_Print(p)
00575 
00576     use globals, only : stream
00577 
00578     type(Polygon) :: p
00579 
00580     integer       :: n, i
00581 
00582 
00583     n = Polygon_getSize(p)
00584     write(stream, '(a,i3,a)') 'Polygon : ',n,' vertices'
00585     write(stream, *) '  sorted : ', Polygon_isSorted(p)
00586     write(stream, *) '  closed : ', Polygon_isClosed(p)
00587 
00588     write(stream, *) '  vertices:'
00589     do i = 1,n
00590        write(stream, '(a,i3,a,e11.5,a,e11.5,a)') '   [',i,'] = (',p%pts%x(i),';',p%pts%y(i),')'
00591     end do
00592   end subroutine Polygon_Print
00593 
00594 end module Polygon_class