grid_service_library/f90/itm__grid__2dpolygon_8f90_source.html

00001 module itm_grid_2dpolygon
00002
00003   use itm_grid_object
00004
00005   implicit none
00006
00007 contains
00008
00009   subroutine get2dPolygonNodes( grid, object, nodes )
00010     type(type_complexgrid), intent(in) :: grid
00011     type(GridObject), intent(in) :: object
00012     type(GridObject), dimension(:), intent(out), allocatable :: nodes
00013
00014     ! internal
00015     type(GridObject), dimension(:), allocatable :: edges, tmpNodes
00016     type(GridObject), dimension(:,:), allocatable :: edgeNodes
00017     integer :: iEdge, iind, ies
00018
00019     ! TODO: Check that object is of the type we expect (i.e., a face)
00020
00021     ! get edges of object
00022     call getComposingObjects( grid, object, edges )
00023
00024     ! get nodes of edges
00025     allocate( edgeNodes(size(edges), 2) )
00026     do iEdge = 1, size(edges)
00027         call getComposingObjects( grid, edges(iEdge), tmpNodes )
00028
00029         !edgeNodes(iEdge, 1:2) = tmpNodes(1:2)
00030         ! Workaround for bug in pgf95
00031         edgeNodes(iEdge, 1) = tmpNodes(1)
00032         edgeNodes(iEdge, 2) = tmpNodes(2)
00033
00034         !edgeNodes(iEdge, 1) = getObject( tmpNodes(1)%cls, tmpNodes(1)%ind )
00035         !edgeNodes(iEdge, 2) = getObject( tmpNodes(2)%cls, tmpNodes(2)%ind )
00036
00037         !call objectAssign( edgeNodes(iEdge, 1), tmpNodes(1) )
00038         !call objectAssign( edgeNodes(iEdge, 2), tmpNodes(2) )
00039     end do
00040
00041     ! sort nodes to form closed boundary of polygon
00042     allocate( nodes(size(edges)+1) )
00043
00044     nodes(1:2) = edgeNodes(1,1:2)
00045     iind = 2
00046
00047     ! Step along the remaining edges
00048     ! First edge already done, don't have to look at it again
00049     do iEdge = 2, size(edges)
00050         ! Find next edge: has node node[iind] and does not have node[iind-1]
00051         do ies = 2, size(edges)
00052             if ( (edgeNodes(ies, 1) == nodes(iind)) &
00053                 & .and. .not. (edgeNodes(ies, 2) == nodes(iind - 1)) ) then
00054                 iind = iind + 1
00055                 nodes(iind) = edgeNodes(ies, 2)
00056                 exit
00057             else if ( (edgeNodes(ies, 2) == nodes(iind)) &
00058                 & .and. .not. (edgeNodes(ies, 1) == nodes(iind - 1))) then
00059                 iind = iind + 1
00060                 nodes(iind) = edgeNodes(ies, 1)
00061                 exit
00062             end if
00063         end do
00064     end do
00065
00066     deallocate(edgeNodes)
00067
00068   end subroutine get2dPolygonNodes
00069
00070
00071   !> Given the coordinates of n nodes defining the boundary
00072   !> of a polygon as argument nodes(n,2),
00073   !> checks whether the point given by argument x(2) is inside the polygon
00074   logical function pointInPolygon( nodes, x )
00075     real(R8), dimension(:,:), intent(in) :: nodes
00076     real(R8), dimension(:), intent(in) :: x
00077
00078     ! internal
00079     integer :: i, cn
00080     real(R8) :: vt
00081
00082     ! check intersections for all nodes
00083     cn = 0
00084     do i = 1, size(nodes, 1) - 1
00085         if (       (nodes(i, 2) <= x(2) .and. nodes(i+1, 2) > x(2)) &
00086             & .or. (nodes(i, 2) > x(2) .and. nodes(i+1, 2) <= x(2)) ) then
00087
00088             vt = (x(2) - nodes(i, 2)) / (nodes(i+1,2) - nodes(i,2))
00089             if ( x(1) < nodes(i,1) + vt * (nodes(i+1,1) - nodes(i,1)) ) then
00090                 cn = cn + 1
00091             end if
00092         end if
00093     end do
00094
00095     ! cn is even: point outside, cn is odd: point inside
00096     pointInPolygon = ( mod(cn,2) == 1 )
00097   end function pointInPolygon
00098
00099
00100 end module itm_grid_2dpolygon