subroutine get_voronoi_area_version_init(g_num,g_xy ,voronoi_area) ! Parameters: ! Input, g_num integer ( kind = 4 ), the number of points. ! g_xy(2,g_num) real ( kind = 8 ) , the coordinates of input points. ! output, voronoi_area(g_num) real ( kind = 8 ) , the area of each voronoi cell . ! the area of the nodes that form the convex hull is flaged by -1.00 ! Local Parameters: ! ! Local, integer ( kind = 4 ) G_DEGREE(G_NUM), the degree of each Voronoi ! cell. ! ! Local, integer ( kind = 4 ) G_FACE(6*G_NUM), the sequence of vertices to ! be used in a traversal of the boundary of the cell associated with each ! point. ! ! Local, integer ( kind = 4 ) G_START(G_NUM), the index in G_FACE of the ! first vertex at which to begin a traversal of the boundary of the ! cell associated with each point. ! ! Local, integer ( kind = 4 ) I_NUM, the number of vertices at infinity of the ! Voronoi diagram. ! ! Local, real ( kind = 8 ) I_XY(2,I_NUM), the direction of the ! vertices at infinity. ! ! Local, integer ( kind = 4 ) V_NUM, the number of vertices of the Voronoi ! diagram. ! ! Local , integer ( kind = 4 ) HULL_NUM, the number of nodes that lie on ! the convex hull. ! ! Local , integer ( kind = 4 ) HULL(G_NUM). Entries 1 through HULL_NUM ! contain the indices of the nodes that form the convex hull, in order. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) g_xy(dim_num,g_num) real ( kind = 8 ) voronoi_area(g_num) integer ( kind = 4 ) g_num integer ( kind = 4 ) hull_num integer ( kind = 4 ) i_num integer ( kind = 4 ) n_num integer ( kind = 4 ) v_num integer ( kind = 4 ) tri_num integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k real ( kind = 8 ) area_temp integer ( kind = 4 ) n_temp integer ( kind = 4 ), allocatable, dimension ( : ) :: g_degree integer ( kind = 4 ), allocatable, dimension ( : ) :: g_face integer ( kind = 4 ), allocatable, dimension ( : ) :: g_start integer ( kind = 4 ), allocatable, dimension ( : ) :: hull integer ( kind = 4 ), allocatable, dimension ( :,: ) :: nodtri integer ( kind = 4 ), allocatable, dimension ( :,: ) :: tnbr integer ( kind = 4 ), allocatable, dimension ( : ) :: indx real ( kind = 8 ), allocatable, dimension ( :, : ) :: v_xy real ( kind = 8 ), allocatable, dimension ( : ) :: cell_x real ( kind = 8 ), allocatable, dimension ( : ) :: cell_y real ( kind = 8 ), allocatable, dimension ( :, : ) :: i_xy allocate ( hull(g_num) ) allocate ( cell_x(20) ) allocate ( cell_y(20) ) allocate ( nodtri(3,2*g_num) ) allocate ( tnbr(3,2*g_num) ) allocate ( indx(g_num) ) call dtris2 ( g_num, g_xy, tri_num, nodtri, tnbr ) call points_hull( g_num,tri_num, nodtri, hull_num, hull ) allocate ( i_xy(dim_num,g_num) ) allocate ( g_degree(g_num) ) allocate ( g_face(10*g_num) ) allocate ( g_start(g_num) ) allocate ( v_xy(dim_num,2*g_num) ) ! get the data defining the Voronoi diagram. call voronoi_data (g_num, g_xy, g_degree, g_start, g_face, v_num, v_xy, & i_num, i_xy ) ! calculate the area of ecah voronoi cell ! open ( unit = 20, file = 'voronoi_xy.txt' ) ! voronoi cell verticals file: do i = 1,g_num k = g_start(i) ! write ( 20, * ) ">> ",i ! for output do j = 1, g_degree(i) cell_x(j) = v_xy(1,g_face(k)) cell_y(j) = v_xy(2,g_face(k)) ! write ( 20, '(2g25.8)' ) v_xy(1,g_face(k)), v_xy(2,g_face(k)) ! for output k = k + 1 end do ! write ( 20,'(2g25.8)') v_xy(1,g_face(g_start(i))),v_xy(2,g_face(g_start(i))) ! for output n_temp = g_degree(i) call areapg ( n_temp , cell_x, cell_y , area_temp ) voronoi_area(i) = area_temp end do ! close ( 20 ,status = "KEEP" ) do i = 1, hull_num voronoi_area(hull(i)) = -1.0 end do deallocate ( i_xy ) deallocate ( g_degree ) deallocate ( g_face ) deallocate ( g_start ) deallocate ( v_xy ) deallocate ( hull ) deallocate ( cell_x ) deallocate ( cell_y ) return end function angle_rad_2d ( p1, p2, p3 ) !*****************************************************************************80 ! !! ANGLE_RAD_2D returns the angle in radians swept out between two rays in 2D. ! ! Discussion: ! ! Except for the zero angle case, it should be true that ! ! ANGLE_RAD_2D ( P1, P2, P3 ) + ANGLE_RAD_2D ( P3, P2, P1 ) = 2 * PI ! ! P1 ! / ! / ! / ! / ! P2--------->P3 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 January 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) P1(2), P2(2), P3(2), define the rays ! P1 - P2 and P3 - P2 which define the angle. ! ! Output, real ( kind = 8 ) ANGLE_RAD_2D, the angle swept out by the rays, ! in radians. 0 <= ANGLE_RAD_2D < 2 * PI. If either ray has zero ! length, then ANGLE_RAD_2D is set to 0. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) angle_rad_2d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) p(dim_num) real ( kind = 8 ) p1(dim_num) real ( kind = 8 ) p2(dim_num) real ( kind = 8 ) p3(dim_num) p(1) = ( p3(1) - p2(1) ) * ( p1(1) - p2(1) ) & + ( p3(2) - p2(2) ) * ( p1(2) - p2(2) ) p(2) = ( p3(1) - p2(1) ) * ( p1(2) - p2(2) ) & - ( p3(2) - p2(2) ) * ( p1(1) - p2(1) ) if ( all ( p(1:dim_num) == 0.0D+00) ) then return end if angle_rad_2d = atan2 ( p(2), p(1) ) if ( angle_rad_2d < 0.0D+00 ) then angle_rad_2d = angle_rad_2d + 2.0D+00 * pi end if return end subroutine angle_half_2d ( p1, p2, p3, p4 ) !*****************************************************************************80 ! !! ANGLE_HALF_2D finds half an angle in 2D. ! ! ! Discussion: ! ! The original angle is defined by the sequence of points P1, P2 and P3. ! ! The point P4 is calculated so that: ! ! (P1,P2,P4) = (P1,P2,P3) / 2 ! ! P1 ! / ! / P4 ! / . ! / . ! P2--------->P3 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 01 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) P1(2), P2(2), P3(2), points defining the angle. ! ! Input, real ( kind = 8 ) P4(2), a point defining the half angle. implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) p1(dim_num) real ( kind = 8 ) p2(dim_num) real ( kind = 8 ) p3(dim_num) real ( kind = 8 ) p4(dim_num) integer ( kind = 4 ) lr integer ( kind = 4 ) lrline lr=lrline(p2(1),p2(2),p1(1),p1(2),p3(1),p3(2),0.0D+00) if ( lr == 0 ) then p4(1)=p2(1)+(p1(2)-p2(2))/(sqrt ( sum ( ( p1(1:2) - p2(1:2) )**2 ) )) p4(2)=p2(2)-(p1(1)-p2(1))/(sqrt ( sum ( ( p1(1:2) - p2(1:2) )**2 ) )) else p4(1:2) = 0.5D+00 * ( & ( p1(1:2) - p2(1:2) ) / sqrt ( sum ( ( p1(1:2) - p2(1:2) )**2 ) ) & + ( p3(1:2) - p2(1:2) ) / sqrt ( sum ( ( p3(1:2) - p2(1:2) )**2 ) ) ) p4(1:2) = p2(1:2) + p4(1:2) / sqrt ( sum ( p4(1:2)**2 ) ) end if return end subroutine voronoi_data ( g_num, g_xy, g_degree, g_start, g_face, v_num, & v_xy, i_num, i_xy ) !*****************************************************************************80 ! !! VORONOI_DATA returns data defining the Voronoi diagram. ! ! Discussion: ! ! The routine first determines the Delaunay triangulation. ! ! The Voronoi diagram is then determined from this information. ! ! In particular, the circumcenter of each Delaunay triangle ! is a vertex of a Voronoi polygon. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) G_NUM, the number of generators. ! ! Input, real ( kind = 8 ) G_XY(2,G_NUM), the point coordinates. ! ! Output, integer ( kind = 4 ) G_DEGREE(G_NUM), the degree of each ! Voronoi cell. ! ! Output, integer ( kind = 4 ) G_START(G_NUM), the index in G_FACE of the ! first vertex at which to begin a traversal of the boundary of the ! cell associated with each point. ! ! Output, integer ( kind = 4 ) G_FACE(6*G_NUM), the sequence of vertices to ! be used in a traversal of the boundary of the cell associated with each ! point. ! ! Output, integer ( kind = 4 ) V_NUM, the number of vertices of the Voronoi ! diagram. ! ! Output, real ( kind = 8 ) V_XY(2,V_NUM), the coordinates of the vertices ! Output, real ( kind = 8 ) I_XY(2,I_NUM), the direction of the ! vertices at infinity. ! implicit none integer ( kind = 4 ) g_num integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) area integer ( kind = 4 ) count logical, parameter :: debug = .true. integer ( kind = 4 ) g integer ( kind = 4 ) g_degree(g_num) integer ( kind = 4 ) g_face(10*g_num) integer ( kind = 4 ) g_next integer ( kind = 4 ) g_start(g_num) real ( kind = 8 ) g_xy(dim_num,g_num) integer ( kind = 4 ) i integer ( kind = 4 ) i_num integer ( kind = 4 ) i4_wrap real ( kind = 8 ) i_xy(dim_num,g_num) integer ( kind = 4 ) i1 integer ( kind = 4 ) i2 integer ( kind = 4 ) i3 integer ( kind = 4 ) ix1 integer ( kind = 4 ) ix2 integer ( kind = 4 ) j integer ( kind = 4 ) jp1 integer ( kind = 4 ) k integer ( kind = 4 ) nodtri(3,2*g_num) integer ( kind = 4 ) s integer ( kind = 4 ) sp1 integer ( kind = 4 ) s_next real ( kind = 8 ) t(dim_num,3) integer ( kind = 4 ) tnbr(3,2*g_num) integer ( kind = 4 ) v integer ( kind = 4 ) v_inf integer ( kind = 4 ) v_next integer ( kind = 4 ) v_num integer ( kind = 4 ) v_old integer ( kind = 4 ) v_save real ( kind = 8 ) v_xy(dim_num,2*g_num) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) y1 real ( kind = 8 ) y2 ! ! Compute the Delaunay triangulation. ! ! do v=1,g_num ! write ( * ,'(2f14.6)' ) g_xy(1,v) , g_xy(2,v) ! end do call dtris2 ( g_num, g_xy, v_num, nodtri, tnbr ) ! ! Compute and print the areas of the finite triangles. ! ! write ( *, '(a)' ) ' ' ! write ( *, '(a)' ) ' Triangle Area' ! write ( *, '(a)' ) ' ' ! do v = 1, v_num i1 = nodtri(1,v) i2 = nodtri(2,v) i3 = nodtri(3,v) t(1:dim_num,1:3) = reshape ( (/ & g_xy(1:2,i1), g_xy(1:2,i2), g_xy(1:2,i3) /), (/ dim_num, 3 /) ) call triangle_area_2d ( t, area ) ! write ( *, '(2x,i8,2x,g14.6)' ) v, area end do ! ! Extend the NODTRI data structure, adding fictitious vertices at infinity, ! so that the Delaunay triangulation can be regarded as covering the ! entire plane. ! call tri_augment ( v_num, nodtri, v_inf ) if ( debug ) then ! write ( *, '(a)' ) ' ' ! write ( *, '(a)' ) ' The generators that form each Delaunay triangle:' ! write ( *, '(a)' ) ' (Negative values are fictitious nodes at infinity.)' ! write ( *, '(a)' ) ' ' ! call i4mat_transpose_print ( 3, v_num+v_inf, nodtri, ' Triangle nodes:' ) end if ! ! Negative entries in TNBR indicate a semi-infinite Voronoi side. ! However, DTRIS2 uses a peculiar numbering. Renumber them. ! i_num = 0 do v = 1, v_num do i = 1, 3 if ( tnbr(i,v) < 0 ) then i_num = i_num + 1 tnbr(i,v) = -i_num end if end do end do ! if ( debug ) then ! write ( *, '(a)' ) ' ' ! write ( *, '(a)' ) ' Neighboring triangles of each Delaunay triangle:' ! write ( *, '(a)' ) ' Negative values indicate no finite neighbor.' ! write ( *, '(a)' ) ' ' ! call i4mat_transpose_print ( 3, v_num, tnbr, ' Neighbor triangles:' ) ! ! end if ! ! Determine the degree of each cell. ! g_degree(1:g_num) = 0 do j = 1, v_num + v_inf do i = 1, 3 k = nodtri(i,j) if ( 0 < k ) then g_degree(k) = g_degree(k) + 1 end if end do end do ! call i4vec_print ( g_num, g_degree, ' Voronoi cell degrees' ) ! ! Each (finite) Delaunay triangle contains a vertex of the Voronoi polygon, ! at the triangle's circumcenter. ! do v = 1, v_num i1 = nodtri(1,v) i2 = nodtri(2,v) i3 = nodtri(3,v) t(1:dim_num,1:3) = reshape ( (/ & g_xy(1:2,i1), g_xy(1:2,i2), g_xy(1:2,i3) /), (/ dim_num, 3 /) ) call triangle_circumcenter_2d ( t, v_xy(1:2,v) ) end do ! call r8mat_transpose_print ( dim_num, v_num, v_xy, ' The Voronoi vertices:' ) ! ! For each generator G: ! Determine if its region is infinite. ! Find a Delaunay triangle containing G. ! Seek another triangle containing the next node in that triangle. ! count = 0 g_start(1:g_num) = 0 do g = 1, g_num v_next = 0 do v = 1, v_num + v_inf do s = 1, 3 if ( nodtri(s,v) == g ) then v_next = v s_next = s exit end if end do if ( v_next /= 0 ) then exit end if end do v_save = v_next do s_next = i4_wrap ( s_next + 1, 1, 3 ) g_next = nodtri(s_next,v_next) if ( g_next == g ) then s_next = i4_wrap ( s_next + 1, 1, 3 ) g_next = nodtri(s_next,v_next) end if v_old = v_next v_next = 0 do v = 1, v_num + v_inf if ( v == v_old ) then cycle end if do s = 1, 3 if ( nodtri(s,v) == g ) then sp1 = i4_wrap ( s + 1, 1, 3 ) if ( nodtri(sp1,v) == g_next ) then v_next = v s_next = sp1 exit end if sp1 = i4_wrap ( s + 2, 1, 3 ) if ( nodtri(sp1,v) == g_next ) then v_next = v s_next = sp1 exit end if end if end do if ( v_next /= 0 ) then exit end if end do if ( v_next == v_save ) then exit end if if ( v_next == 0 ) then v_next = v_old exit end if end do ! ! Now, starting in the current triangle, V_NEXT, cycle again, ! and copy the list of nodes into the array. ! v_save = v_next count = count + 1 g_start(g) = count g_face(count) = v_next do s_next = i4_wrap ( s_next + 1, 1, 3 ) g_next = nodtri(s_next,v_next) if ( g_next == g ) then s_next = i4_wrap ( s_next + 1, 1, 3 ) g_next = nodtri(s_next,v_next) end if v_old = v_next v_next = 0 do v = 1, v_num + v_inf if ( v == v_old ) then cycle end if do s = 1, 3 if ( nodtri(s,v) == g ) then sp1 = i4_wrap ( s + 1, 1, 3 ) if ( nodtri(sp1,v) == g_next ) then v_next = v s_next = sp1 exit end if sp1 = i4_wrap ( s + 2, 1, 3 ) if ( nodtri(sp1,v) == g_next ) then v_next = v s_next = sp1 exit end if end if end do if ( v_next /= 0 ) then exit end if end do if ( v_next == v_save ) then exit end if if ( v_next == 0 ) then exit end if count = count + 1 g_face(count) = v_next end do end do ! ! Mark all the vertices at infinity with a negative sign, ! so that the data in G_FACE is easier to interpret. ! do i = 1, count if ( v_num < g_face(i) ) then g_face(i) = -g_face(i) end if end do ! ! For each (finite) Delaunay triangle, I ! For each side J, ! do i = 1, v_num do j = 1, 3 k = tnbr(j,i) ! If there is no neighboring triangle on that side, ! extend a line from the circumcenter of I in the direction of the ! outward normal to that side. This is an infinite edge of ! an infinite Voronoi polygon. ! if ( k < 0 ) then ix1 = nodtri(j,i) x1 = g_xy(1,ix1) y1 = g_xy(2,ix1) jp1 = i4_wrap ( j+1, 1, 3 ) ix2 = nodtri(jp1,i) x2 = g_xy(1,ix2) y2 = g_xy(2,ix2) ! ! Compute the direction I_XY(1:2,-K). ! ! call line_exp_normal_2d ( g_xy(1:2,ix1), g_xy(1:2,ix2), i_xy(1:2,-k) ) end if end do end do return end subroutine dtris2 ( point_num, point_xy, tri_num, tri_vert, tri_nabe ) !*****************************************************************************80 ! !! DTRIS2 constructs a Delaunay triangulation of 2D vertices. ! ! Discussion: ! ! The routine constructs the Delaunay triangulation of a set of 2D vertices ! using an incremental approach and diagonal edge swaps. Vertices are ! first sorted in lexicographically increasing (X,Y) order, and ! then are inserted one at a time from outside the convex hull. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 January 2004 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of vertices. ! ! Input/output, real ( kind = 8 ) POINT_XY(2,POINT_NUM), the vertices. ! On output, the vertices have been sorted into dictionary order. ! ! Output, integer ( kind = 4 ) TRI_NUM, the number of triangles in the ! triangulation; TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the ! number of boundary vertices. ! ! Output, integer ( kind = 4 ) TRI_VERT(3,TRI_NUM), the nodes that make up ! each triangle. The elements are indices of POINT_XY. The vertices of the ! triangles are in counter clockwise order. ! ! ! Output, integer ( kind = 4 ) TRI_NABE(3,TRI_NUM), the triangle neighbor ! list. Positive elements are indices of TIL; negative elements are used ! for links of a counter clockwise linked list of boundary edges; ! LINK = -(3*I + J-1) where I, J = triangle, edge index; TRI_NABE(J,I) refers ! to the neighbor along edge from vertex J to J+1 (mod 3). ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) point_num real ( kind = 8 ) cmax integer ( kind = 4 ) e integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) indx(point_num) integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) ledg integer ( kind = 4 ) lr integer ( kind = 4 ) lrline integer ( kind = 4 ) ltri integer ( kind = 4 ) m integer ( kind = 4 ) m1 integer ( kind = 4 ) m2 integer ( kind = 4 ) n real ( kind = 8 ) point_xy(dim_num,point_num) integer ( kind = 4 ) redg integer ( kind = 4 ) rtri integer ( kind = 4 ) stack(point_num) integer ( kind = 4 ) t real ( kind = 8 ) tol integer ( kind = 4 ) top integer ( kind = 4 ) tri_nabe(3,point_num*2) integer ( kind = 4 ) tri_num integer ( kind = 4 ) tri_vert(3,point_num*2) tol = 100.0D+00 * epsilon ( tol ) ierr = 0 ! ! Sort the vertices by increasing (x,y). ! call r82vec_sort_heap_index_a ( point_num, point_xy, indx ) call r82vec_permute ( point_num, point_xy, indx ) ! ! Make sure that the data points are "reasonably" distinct. ! m1 = 1 do i = 2, point_num m = m1 m1 = i k = 0 do j = 1, 2 cmax = max ( abs ( point_xy(j,m) ), abs ( point_xy(j,m1) ) ) if ( tol * ( cmax + 1.0D+00 ) & < abs ( point_xy(j,m) - point_xy(j,m1) ) ) then k = j exit end if end do if ( k == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a,i8)' ) ' Fails for point number I = ', i write ( *, '(a,i8)' ) ' M = ', m write ( *, '(a,i8)' ) ' M1 = ', m1 write ( *, '(a,2g14.6)' ) ' X,Y(M) = ', point_xy(1,m), point_xy(2,m) write ( *, '(a,2g14.6)' ) ' X,Y(M1) = ', point_xy(1,m1), point_xy(2,m1) ierr = 224 return end if end do ! ! Starting from points M1 and M2, search for a third point M that ! makes a "healthy" triangle (M1,M2,M) ! m1 = 1 m2 = 2 j = 3 do if ( point_num < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' ierr = 225 return end if m = j lr = lrline ( point_xy(1,m), point_xy(2,m), point_xy(1,m1), & point_xy(2,m1), point_xy(1,m2), point_xy(2,m2), 0.0D+00 ) if ( lr /= 0 ) then exit end if j = j + 1 end do ! ! Set up the triangle information for (M1,M2,M), and for any other ! triangles you created because points were collinear with M1, M2. ! tri_num = j - 2 if ( lr == -1 ) then tri_vert(1,1) = m1 tri_vert(2,1) = m2 tri_vert(3,1) = m tri_nabe(3,1) = -3 do i = 2, tri_num m1 = m2 m2 = i+1 tri_vert(1,i) = m1 tri_vert(2,i) = m2 tri_vert(3,i) = m tri_nabe(1,i-1) = -3 * i tri_nabe(2,i-1) = i tri_nabe(3,i) = i - 1 end do tri_nabe(1,tri_num) = -3 * tri_num - 1 tri_nabe(2,tri_num) = -5 ledg = 2 ltri = tri_num else tri_vert(1,1) = m2 tri_vert(2,1) = m1 tri_vert(3,1) = m tri_nabe(1,1) = -4 do i = 2, tri_num m1 = m2 m2 = i+1 tri_vert(1,i) = m2 tri_vert(2,i) = m1 tri_vert(3,i) = m tri_nabe(3,i-1) = i tri_nabe(1,i) = -3 * i - 3 tri_nabe(2,i) = i - 1 end do tri_nabe(3,tri_num) = -3 * tri_num tri_nabe(2,1) = -3 * tri_num - 2 ledg = 2 ltri = 1 end if ! ! ! Insert the vertices one at a time from outside the convex hull, ! determine visible boundary edges, and apply diagonal edge swaps until ! Delaunay triangulation of vertices (so far) is obtained. ! top = 0 do i = j+1, point_num m = i m1 = tri_vert(ledg,ltri) if ( ledg <= 2 ) then m2 = tri_vert(ledg+1,ltri) else m2 = tri_vert(1,ltri) end if lr = lrline ( point_xy(1,m), point_xy(2,m), point_xy(1,m1), & point_xy(2,m1), point_xy(1,m2), point_xy(2,m2), 0.0D+00 ) if ( 0 < lr ) then rtri = ltri redg = ledg ltri = 0 else l = -tri_nabe(ledg,ltri) rtri = l / 3 redg = mod(l,3) + 1 end if call vbedg ( point_xy(1,m), point_xy(2,m), point_num, point_xy, tri_num, & tri_vert, tri_nabe, ltri, ledg, rtri, redg ) n = tri_num + 1 l = -tri_nabe(ledg,ltri) do t = l / 3 e = mod ( l, 3 ) + 1 l = -tri_nabe(e,t) m2 = tri_vert(e,t) if ( e <= 2 ) then m1 = tri_vert(e+1,t) else m1 = tri_vert(1,t) end if tri_num = tri_num + 1 tri_nabe(e,t) = tri_num tri_vert(1,tri_num) = m1 tri_vert(2,tri_num) = m2 tri_vert(3,tri_num) = m tri_nabe(1,tri_num) = t tri_nabe(2,tri_num) = tri_num - 1 tri_nabe(3,tri_num) = tri_num + 1 top = top + 1 if ( point_num < top ) then ierr = 8 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a)' ) ' Stack overflow.' return end if stack(top) = tri_num if ( t == rtri .and. e == redg ) then exit end if end do tri_nabe(ledg,ltri) = -3 * n - 1 tri_nabe(2,n) = -3 * tri_num - 2 tri_nabe(3,tri_num) = -l ltri = n ledg = 2 call swapec ( m, top, ltri, ledg, point_num, point_xy, tri_num, & tri_vert, tri_nabe, stack, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a)' ) ' Error return from SWAPEC.' return end if end do ! ! Now account for the sorting that we did. ! do i = 1, 3 do j = 1, tri_num tri_vert(i,j) = indx ( tri_vert(i,j) ) end do end do call perm_inv ( point_num, indx ) call r82vec_permute ( point_num, point_xy, indx ) return end subroutine areapg ( nvrt, xc, yc, area ) implicit none integer ( kind = 4 ) nvrt real ( kind = 8 ) area integer ( kind = 4 ) i real ( kind = 8 ) sum2 real ( kind = 8 ) xc(10) real ( kind = 8 ) yc(10) sum2 = xc(1) * ( yc(2) - yc(nvrt) ) do i = 2, nvrt-1 sum2 = sum2 + xc(i) * ( yc(i+1) - yc(i-1) ) end do sum2 = sum2 + xc(nvrt) * ( yc(1) - yc(nvrt-1) ) area = ABS(sum2*0.5) ! write( *,'(a)' ) 'area cal sus~' return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an I4 to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I4_WRAP ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) IVAL, an integer value. ! ! Input, integer ( kind = 4 ) ILO, IHI, the desired bounds. ! ! Output, integer ( kind = 4 ) I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer ( kind = 4 ) i4_modp integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) ival integer ( kind = 4 ) wide wide = ihi + 1 - ilo if ( wide == 0 ) then i4_wrap = ilo else i4_wrap = ilo + i4_modp ( ival - ilo, wide ) end if return end function lrline ( xu, yu, xv1, yv1, xv2, yv2, dv ) !*****************************************************************************80 ! !! LRLINE determines if a point is left of, right or, or on a directed line. ! ! Discussion: ! ! The directed line is parallel to, and at a signed distance DV from ! a directed base line from (XV1,YV1) to (XV2,YV2). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) XU, YU, the coordinates of the point whose ! position relative to the directed line is to be determined. ! ! Input, real ( kind = 8 ) XV1, YV1, XV2, YV2, the coordinates of two points ! that determine the directed base line. ! ! Input, real ( kind = 8 ) DV, the signed distance of the directed line ! from the directed base line through the points (XV1,YV1) and (XV2,YV2). ! DV is positive for a line to the left of the base line. ! ! Output, integer ( kind = 4 ) LRLINE, the result: ! +1, the point is to the right of the directed line; ! 0, the point is on the directed line; ! -1, the point is to the left of the directed line. ! implicit none real ( kind = 8 ) dv real ( kind = 8 ) dx real ( kind = 8 ) dxu real ( kind = 8 ) dy real ( kind = 8 ) dyu integer ( kind = 4 ) lrline real ( kind = 8 ) t real ( kind = 8 ) tol real ( kind = 8 ) tolabs real ( kind = 8 ) xu real ( kind = 8 ) xv1 real ( kind = 8 ) xv2 real ( kind = 8 ) yu real ( kind = 8 ) yv1 real ( kind = 8 ) yv2 tol = 100.0D+00 * epsilon ( tol ) dx = xv2 - xv1 dy = yv2 - yv1 dxu = xu - xv1 dyu = yu - yv1 tolabs = tol * max ( abs ( dx ), abs ( dy ), abs ( dxu ), & abs ( dyu ), abs ( dv ) ) t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy ) if ( tolabs < t ) then lrline = 1 else if ( -tolabs <= t ) then lrline = 0 else lrline = -1 end if return end subroutine tri_augment ( v_num, nodtri, v_inf ) !*****************************************************************************80 ! !! TRI_AUGMENT augments the triangle data using vertices at infinity. ! ! Discussion: ! ! The algorithm simply looks at the list of triangle edges stored ! in NODTRI, and determines which edges, of the form (P1,P2), do ! not have a matching (P2,P1) occurrence. These correspond to ! boundary edges of the convex hull. To simplify our computations, ! we adjust the NODTRI array to accommodate an extra triangle with ! one vertex at infinity for each such unmatched edge. ! ! The algorithm used here is ruinously inefficient for large V_NUM. ! Assuming that this data structure modification is the way to go, ! the routine should be rewritten to determine the boundary edges ! more efficiently. ! ! The fictitious vertices at infinity show up in the augmenting ! rows of the NODTRI array with negative indices. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) V_NUM, the number of Voronoi vertices. ! ! Input/output, integer ( kind = 4 ) NODTRI(3,*), the list of nodes that ! comprise each Delaunay triangle. On input, there are V_NUM ! sets of this data. On output, for every pair of nodes (P1,P2) ! for which the pair (P2,P1) does not occur, an augmenting triangle ! has been created with exactly this edge (plus a vertex at infinity). ! On output, there are V_NUM + V_INF sets of data. ! ! Output, integer ( kind = 4 ) V_INF, the number of augmenting triangles and ! vertices at infinity that were created. ! implicit none logical found integer ( kind = 4 ) i integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) i2 integer ( kind = 4 ) ip1 integer ( kind = 4 ) nodtri(3,*) integer ( kind = 4 ) s integer ( kind = 4 ) s2 integer ( kind = 4 ) t integer ( kind = 4 ) t2 integer ( kind = 4 ) v integer ( kind = 4 ) v_inf integer ( kind = 4 ) v_num integer ( kind = 4 ) v2 v_inf = 0 do v = 1, v_num do i = 1, 3 s = nodtri(i,v) ip1 = i4_wrap ( i + 1, 1, 3 ) t = nodtri(ip1,v) found = .false. do v2 = 1, v_num do i2 = 1, 3 s2 = nodtri(i2,v2) ip1 = i4_wrap ( i2 + 1, 1, 3 ) t2 = nodtri(ip1,v2) if ( s == t2 .and. t == s2 ) then found = .true. exit end if end do if ( found ) then exit end if end do if ( .not. found ) then v_inf = v_inf + 1 nodtri(1:3,v_num+v_inf) = (/ -v_inf, t, s /) end if end do end do ! write ( *, '(a)' ) ' ' ! write ( *, '(a)' ) 'TRI_AUGMENT:' ! write ( *, '(a,i8)' ) ' Number of boundary triangles = ', v_inf return end subroutine triangle_circumcenter_2d ( t, center ) !*****************************************************************************80 ! !! TRIANGLE_CIRCUMCENTER_2D computes the circumcenter of a triangle in 2D. ! ! Discussion: ! ! The circumcenter of a triangle is the center of the circumcircle, the ! circle that passes through the three vertices of the triangle. ! ! The circumcircle contains the triangle, but it is not necessarily the ! smallest triangle to do so. ! ! If all angles of the triangle are no greater than 90 degrees, then ! the center of the circumscribed circle will lie inside the triangle. ! Otherwise, the center will lie outside the triangle. ! ! The circumcenter is the intersection of the perpendicular bisectors ! of the sides of the triangle. ! ! In geometry, the circumcenter of a triangle is often symbolized by "O". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,3), the triangle vertices. ! ! Output, real ( kind = 8 ) CENTER(2), the circumcenter of the triangle. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) asq real ( kind = 8 ) bot real ( kind = 8 ) center(dim_num) real ( kind = 8 ) csq real ( kind = 8 ) t(dim_num,3) real ( kind = 8 ) top(dim_num) asq = ( t(1,2) - t(1,1) )**2 + ( t(2,2) - t(2,1) )**2 csq = ( t(1,3) - t(1,1) )**2 + ( t(2,3) - t(2,1) )**2 top(1) = ( t(2,2) - t(2,1) ) * csq - ( t(2,3) - t(2,1) ) * asq top(2) = - ( t(1,2) - t(1,1) ) * csq + ( t(1,3) - t(1,1) ) * asq bot = ( t(2,2) - t(2,1) ) * ( t(1,3) - t(1,1) ) & - ( t(2,3) - t(2,1) ) * ( t(1,2) - t(1,1) ) center(1:2) = t(1:2,1) + 0.5D+00 * top(1:2) / bot return end subroutine vbedg ( x, y, point_num, point_xy, tri_num, tri_vert, tri_nabe, & ltri, ledg, rtri, redg ) implicit none integer ( kind = 4 ) point_num integer ( kind = 4 ) tri_num integer ( kind = 4 ) a integer ( kind = 4 ) b integer ( kind = 4 ) e integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) l logical ldone integer ( kind = 4 ) ledg integer ( kind = 4 ) lr integer ( kind = 4 ) lrline integer ( kind = 4 ) ltri real ( kind = 8 ) point_xy(2,point_num) integer ( kind = 4 ) redg integer ( kind = 4 ) rtri integer ( kind = 4 ) t integer ( kind = 4 ) tri_nabe(3,tri_num) integer ( kind = 4 ) tri_vert(3,tri_num) real ( kind = 8 ) x real ( kind = 8 ) y ! ! Find the rightmost visible boundary edge using links, then possibly ! leftmost visible boundary edge using triangle neighbor information. ! if ( ltri == 0 ) then ldone = .false. ltri = rtri ledg = redg else ldone = .true. end if do l = -tri_nabe(redg,rtri) t = l / 3 e = mod ( l, 3 ) + 1 a = tri_vert(e,t) if ( e <= 2 ) then b = tri_vert(e+1,t) else b = tri_vert(1,t) end if lr = lrline ( x, y, point_xy(1,a), point_xy(2,a), point_xy(1,b), & point_xy(2,b), 0.0D+00 ) if ( lr <= 0 ) then exit end if rtri = t redg = e end do if ( ldone ) then return end if t = ltri e = ledg do b = tri_vert(e,t) e = i4_wrap ( e-1, 1, 3 ) do while ( 0 < tri_nabe(e,t) ) t = tri_nabe(e,t) if ( tri_vert(1,t) == b ) then e = 3 else if ( tri_vert(2,t) == b ) then e = 1 else e = 2 end if end do a = tri_vert(e,t) lr = lrline ( x, y, point_xy(1,a), point_xy(2,a), point_xy(1,b), & point_xy(2,b), 0.0D+00 ) if ( lr <= 0 ) then exit end if end do ltri = t ledg = e return end function i4_modp ( i, j ) implicit none integer ( kind = 4 ) i integer ( kind = 4 ) i4_modp integer ( kind = 4 ) j if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i8)' ) ' I4_MODP ( I, J ) called with J = ', j stop end if i4_modp = mod ( i, j ) if ( i4_modp < 0 ) then i4_modp = i4_modp + abs ( j ) end if return end subroutine i4vec_indicator ( n, a ) !*****************************************************************************80 ! !! I4VEC_INDICATOR sets an I4VEC to the indicator vector. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of elements of A. ! ! Output, integer ( kind = 4 ) A(N), the array to be initialized. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) i do i = 1, n a(i) = i end do return end subroutine triangle_area_2d ( t, area ) !*****************************************************************************80 ! !! TRIANGLE_AREA_2D computes the area of a triangle in 2D. ! ! Discussion: ! ! If the triangle's vertices are given in counter clockwise order, ! the area will be positive. If the triangle's vertices are given ! in clockwise order, the area will be negative! ! ! An earlier version of this routine always returned the absolute ! value of the computed area. I am convinced now that that is ! a less useful result! For instance, by returning the signed ! area of a triangle, it is possible to easily compute the area ! of a nonconvex polygon as the sum of the (possibly negative) ! areas of triangles formed by node 1 and successive pairs of vertices. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,3), the triangle vertices. ! ! Output, real ( kind = 8 ) AREA, the area of the triangle. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) area real ( kind = 8 ) t(dim_num,3) area = 0.5D+00 * ( & t(1,1) * ( t(2,2) - t(2,3) ) & + t(1,2) * ( t(2,3) - t(2,1) ) & + t(1,3) * ( t(2,1) - t(2,2) ) ) return end subroutine points_hull( g_num,v_num, nodtri, v_inf, hull ) implicit none logical found integer ( kind = 4 ) i integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) i2 integer ( kind = 4 ) ip1 integer ( kind = 4 ) nodtri(3,*) integer ( kind = 4 ) s integer ( kind = 4 ) s2 integer ( kind = 4 ) t integer ( kind = 4 ) t2 integer ( kind = 4 ) v integer ( kind = 4 ) v_inf integer ( kind = 4 ) v_num integer ( kind = 4 ) v2 integer ( kind = 4 ) g_num integer ( kind = 4 ) hull(g_num) logical point_flag(g_num) integer ( kind = 4 ) k v_inf = 0 do i = 1,g_num point_flag(i)=.false. end do do v = 1, v_num do i = 1, 3 s = nodtri(i,v) ip1 = i4_wrap ( i + 1, 1, 3 ) t = nodtri(ip1,v) found = .false. do v2 = 1, v_num do i2 = 1, 3 s2 = nodtri(i2,v2) ip1 = i4_wrap ( i2 + 1, 1, 3 ) t2 = nodtri(ip1,v2) if ( s == t2 .and. t == s2 ) then found = .true. exit end if end do if ( found ) then exit end if end do if ( .not. found ) then v_inf = v_inf + 1 ! nodtri(1:3,v_num+v_inf) = (/ -v_inf, t, s /) point_flag(t) = .true. point_flag(s) = .true. end if end do end do k = 0 do i = 1,g_num if ( point_flag(i) ) then k= k+1 hull(k) = i end if end do return end subroutine r82vec_permute ( n, a, p ) !*****************************************************************************80 ! !! R82VEC_PERMUTE permutes an R82VEC in place. ! ! Discussion: ! ! This routine permutes an array of real "objects", but the same ! logic can be used to permute an array of objects of any arithmetic ! type, or an array of objects of any complexity. The only temporary ! storage required is enough to store a single object. The number ! of data movements made is N + the number of cycles of order 2 or more, ! which is never more than N + N/2. ! ! Example: ! ! Input: ! ! N = 5 ! P = ( 2, 4, 5, 1, 3 ) ! A = ( 1.0, 2.0, 3.0, 4.0, 5.0 ) ! (11.0, 22.0, 33.0, 44.0, 55.0 ) ! ! Output: ! ! A = ( 2.0, 4.0, 5.0, 1.0, 3.0 ) ! ( 22.0, 44.0, 55.0, 11.0, 33.0 ). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of objects. ! ! Input/output, real ( kind = 8 ) A(2,N), the array to be permuted. ! ! Input, integer ( kind = 4 ) P(N), the permutation. P(I) = J means ! that the I-th element of the output array should be the J-th ! element of the input array. P must be a legal permutation ! of the integers from 1 to N, otherwise the algorithm will ! fail catastrophically. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(2,n) real ( kind = 8 ) a_temp(2) integer ( kind = 4 ) iget integer ( kind = 4 ) iput integer ( kind = 4 ) istart integer ( kind = 4 ) p(n) ! ! Search for the next element of the permutation that has not been used. ! do istart = 1, n if ( p(istart) < 0 ) then cycle else if ( p(istart) == istart ) then p(istart) = - p(istart) cycle else a_temp(1:2) = a(1:2,istart) iget = istart ! ! Copy the new value into the vacated entry. ! do iput = iget iget = p(iget) p(iput) = - p(iput) if ( iget < 1 .or. n < iget ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R82VEC_PERMUTE - Fatal error!' stop end if if ( iget == istart ) then a(1:2,iput) = a_temp(1:2) exit end if a(1:2,iput) = a(1:2,iget) end do end if end do ! ! Restore the signs of the entries. ! p(1:n) = -p(1:n) return end subroutine r82vec_sort_heap_index_a ( n, a, indx ) !*****************************************************************************80 ! !! R82VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R82VEC. ! ! Discussion: ! ! The sorting is not actually carried out. Rather an index array is ! created which defines the sorting. This array may be used to sort ! or index the array, or to sort or index related arrays keyed on the ! original array. ! ! Once the index array is computed, the sorting can be carried out ! "implicitly: ! ! A(1:2,INDX(I)), I = 1 to N is sorted, ! ! or explicitly, by the call ! ! call R82VEC_PERMUTE ( N, A, INDX ) ! ! after which A(1:2,I), I = 1 to N is sorted. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the array. ! ! Input, real ( kind = 8 ) A(2,N), an array to be index-sorted. ! ! Output, integer ( kind = 4 ) INDX(N), the sort index. The ! I-th element of the sorted array is A(1:2,INDX(I)). ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(2,n) real ( kind = 8 ) aval(2) integer ( kind = 4 ) i integer ( kind = 4 ) indx(n) integer ( kind = 4 ) indxt integer ( kind = 4 ) ir integer ( kind = 4 ) j integer ( kind = 4 ) l if ( n < 1 ) then return end if if ( n == 1 ) then indx(1) = 1 return end if call i4vec_indicator ( n, indx ) l = n / 2 + 1 ir = n do if ( 1 < l ) then l = l - 1 indxt = indx(l) aval(1:2) = a(1:2,indxt) else indxt = indx(ir) aval(1:2) = a(1:2,indxt) indx(ir) = indx(1) ir = ir - 1 if ( ir == 1 ) then indx(1) = indxt exit end if end if i = l j = l + l do while ( j <= ir ) if ( j < ir ) then if ( a(1,indx(j)) < a(1,indx(j+1)) .or. & ( a(1,indx(j)) == a(1,indx(j+1)) .and. & a(2,indx(j)) < a(2,indx(j+1)) ) ) then j = j + 1 end if end if if ( aval(1) < a(1,indx(j)) .or. & ( aval(1) == a(1,indx(j)) .and. & aval(2) < a(2,indx(j)) ) ) then indx(i) = indx(j) i = j j = j + j else j = ir + 1 end if end do indx(i) = indxt end do return end subroutine swapec ( i, top, btri, bedg, point_num, point_xy, tri_num, & tri_vert, tri_nabe, stack, ierr ) !*****************************************************************************80 ! !! SWAPEC swaps diagonal edges until all triangles are Delaunay. ! ! Discussion: ! ! The routine swaps diagonal edges in a 2D triangulation, based on ! the empty circumcircle criterion, until all triangles are Delaunay, ! given that I is the index of the new vertex added to the triangulation. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) I, the index of the new vertex. ! ! Input/output, integer ( kind = 4 ) TOP, the index of the top of the stack. ! On output, TOP is zero. ! ! Input/output, integer ( kind = 4 ) BTRI, BEDG; on input, if positive, are ! the triangle and edge indices of a boundary edge whose updated indices ! must be recorded. On output, these may be updated because of swaps. ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of points. ! ! Input, real ( kind = 8 ) POINT_XY(2,POINT_NUM), the coordinates ! of the points. ! ! Input, integer ( kind = 4 ) TRI_NUM, the number of triangles. ! ! Input/output, integer ( kind = 4 ) TRI_VERT(3,TRI_NUM), the triangle ! incidence list. May be updated on output because of swaps. ! ! Input/output, integer ( kind = 4 ) TRI_NABE(3,TRI_NUM), the triangle ! neighbor list; negative values are used for links of the counter-clockwise ! linked list of boundary edges; May be updated on output because of swaps. ! ! LINK = -(3*I + J-1) where I, J = triangle, edge index. ! ! Workspace, integer ( kind = 4 ) STACK(MAXST); on input, entries 1 through ! TOP contain the indices of initial triangles (involving vertex I) ! put in stack; the edges opposite I should be in interior; entries ! TOP+1 through MAXST are used as a stack. ! ! Output, integer IERR is set to 8 for abnormal return. ! implicit none integer ( kind = 4 ) point_num integer ( kind = 4 ) tri_num integer ( kind = 4 ) a integer ( kind = 4 ) b integer ( kind = 4 ) bedg integer ( kind = 4 ) btri integer ( kind = 4 ) c integer ( kind = 4 ) diaedg integer ( kind = 4 ) e integer ( kind = 4 ) ee integer ( kind = 4 ) em1 integer ( kind = 4 ) ep1 integer ( kind = 4 ) f integer ( kind = 4 ) fm1 integer ( kind = 4 ) fp1 integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) l real ( kind = 8 ) point_xy(2,point_num) integer ( kind = 4 ) r integer ( kind = 4 ) s integer ( kind = 4 ) stack(point_num) integer ( kind = 4 ) swap integer ( kind = 4 ) t integer ( kind = 4 ) top integer ( kind = 4 ) tri_nabe(3,tri_num) integer ( kind = 4 ) tri_vert(3,tri_num) integer ( kind = 4 ) tt integer ( kind = 4 ) u real ( kind = 8 ) x real ( kind = 8 ) y ! ! Determine whether triangles in stack are Delaunay, and swap ! diagonal edge of convex quadrilateral if not. ! x = point_xy(1,i) y = point_xy(2,i) do if ( top <= 0 ) then exit end if t = stack(top) top = top - 1 if ( tri_vert(1,t) == i ) then e = 2 b = tri_vert(3,t) else if ( tri_vert(2,t) == i ) then e = 3 b = tri_vert(1,t) else e = 1 b = tri_vert(2,t) end if a = tri_vert(e,t) u = tri_nabe(e,t) if ( tri_nabe(1,u) == t ) then f = 1 c = tri_vert(3,u) else if ( tri_nabe(2,u) == t ) then f = 2 c = tri_vert(1,u) else f = 3 c = tri_vert(2,u) end if swap = diaedg ( x, y, point_xy(1,a), point_xy(2,a), point_xy(1,c), & point_xy(2,c), point_xy(1,b), point_xy(2,b) ) if ( swap == 1 ) then em1 = i4_wrap ( e - 1, 1, 3 ) ep1 = i4_wrap ( e + 1, 1, 3 ) fm1 = i4_wrap ( f - 1, 1, 3 ) fp1 = i4_wrap ( f + 1, 1, 3 ) tri_vert(ep1,t) = c tri_vert(fp1,u) = i r = tri_nabe(ep1,t) s = tri_nabe(fp1,u) tri_nabe(ep1,t) = u tri_nabe(fp1,u) = t tri_nabe(e,t) = s tri_nabe(f,u) = r if ( 0 < tri_nabe(fm1,u) ) then top = top + 1 stack(top) = u end if if ( 0 < s ) then if ( tri_nabe(1,s) == u ) then tri_nabe(1,s) = t else if ( tri_nabe(2,s) == u ) then tri_nabe(2,s) = t else tri_nabe(3,s) = t end if top = top + 1 if ( point_num < top ) then ierr = 8 return end if stack(top) = t else if ( u == btri .and. fp1 == bedg ) then btri = t bedg = e end if l = - ( 3 * t + e - 1 ) tt = t ee = em1 do while ( 0 < tri_nabe(ee,tt) ) tt = tri_nabe(ee,tt) if ( tri_vert(1,tt) == a ) then ee = 3 else if ( tri_vert(2,tt) == a ) then ee = 1 else ee = 2 end if end do tri_nabe(ee,tt) = l end if if ( 0 < r ) then if ( tri_nabe(1,r) == t ) then tri_nabe(1,r) = u else if ( tri_nabe(2,r) == t ) then tri_nabe(2,r) = u else tri_nabe(3,r) = u end if else if ( t == btri .and. ep1 == bedg ) then btri = u bedg = f end if l = - ( 3 * u + f - 1 ) tt = u ee = fm1 do while ( 0 < tri_nabe(ee,tt) ) tt = tri_nabe(ee,tt) if ( tri_vert(1,tt) == b ) then ee = 3 else if ( tri_vert(2,tt) == b ) then ee = 1 else ee = 2 end if end do tri_nabe(ee,tt) = l end if end if end do return end function diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 ) !*****************************************************************************80 ! !! DIAEDG chooses a diagonal edge. ! ! Discussion: ! ! The routine determines whether 0--2 or 1--3 is the diagonal edge ! that should be chosen, based on the circumcircle criterion, where ! (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple ! quadrilateral in counterclockwise order. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) X0, Y0, X1, Y1, X2, Y2, X3, Y3, the ! coordinates of the vertices of a quadrilateral, given in ! counter clockwise order. ! ! Output, integer ( kind = 4 ) DIAEDG, chooses a diagonal: ! +1, if diagonal edge 02 is chosen; ! -1, if diagonal edge 13 is chosen; ! 0, if the four vertices are cocircular. ! implicit none real ( kind = 8 ) ca real ( kind = 8 ) cb integer ( kind = 4 ) diaedg real ( kind = 8 ) dx10 real ( kind = 8 ) dx12 real ( kind = 8 ) dx30 real ( kind = 8 ) dx32 real ( kind = 8 ) dy10 real ( kind = 8 ) dy12 real ( kind = 8 ) dy30 real ( kind = 8 ) dy32 real ( kind = 8 ) s real ( kind = 8 ) tol real ( kind = 8 ) tola real ( kind = 8 ) tolb real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y0 real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 tol = 100.0D+00 * epsilon ( tol ) dx10 = x1 - x0 dy10 = y1 - y0 dx12 = x1 - x2 dy12 = y1 - y2 dx30 = x3 - x0 dy30 = y3 - y0 dx32 = x3 - x2 dy32 = y3 - y2 tola = tol * max ( abs ( dx10 ), abs ( dy10 ), abs ( dx30 ), abs ( dy30 ) ) tolb = tol * max ( abs ( dx12 ), abs ( dy12 ), abs ( dx32 ), abs ( dy32 ) ) ca = dx10 * dx30 + dy10 * dy30 cb = dx12 * dx32 + dy12 * dy32 if ( tola < ca .and. tolb < cb ) then diaedg = -1 else if ( ca < -tola .and. cb < -tolb ) then diaedg = 1 else tola = max ( tola, tolb ) s = ( dx10 * dy30 - dx30 * dy10 ) * cb + ( dx32 * dy12 - dx12 * dy32 ) * ca if ( tola < s ) then diaedg = -1 else if ( s < -tola ) then diaedg = 1 else diaedg = 0 end if end if return end subroutine perm_inv ( n, p ) !*****************************************************************************80 ! !! PERM_INV inverts a permutation "in place". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of objects being permuted. ! ! Input/output, integer ( kind = 4 ) P(N), the permutation, in standard ! index form. On output, P describes the inverse permutation ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) i0 integer ( kind = 4 ) i1 integer ( kind = 4 ) i2 integer ( kind = 4 ) is integer ( kind = 4 ) p(n) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PERM_INV - Fatal error!' write ( *, '(a,i8)' ) ' Input value of N = ', n stop end if is = 1 do i = 1, n i1 = p(i) do while ( i < i1 ) i2 = p(i1) p(i1) = -i2 i1 = i2 end do is = -sign ( 1, p(i) ) p(i) = sign ( p(i), is ) end do do i = 1, n i1 = -p(i) if ( 0 <= i1 ) then i0 = i do i2 = p(i1) p(i1) = i0 if ( i2 < 0 ) then exit end if i0 = i1 i1 = i2 end do end if end do return end subroutine ch_cap ( c ) implicit none character c integer ( kind = 4 ) itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end