function alpha_measure ( n, z, triangle_order, triangle_num, triangle_node ) !*****************************************************************************80 ! !! ALPHA_MEASURE determines the triangulated pointset quality measure ALPHA. ! ! Discusion: ! ! The ALPHA measure evaluates the uniformity of the shapes of the triangles ! defined by a triangulated pointset. ! ! We compute the minimum angle among all the triangles in the triangulated ! dataset and divide by the maximum possible value (which, in degrees, ! is 60). The best possible value is 1, and the worst 0. A good ! triangulation should have an ALPHA score close to 1. ! ! The code has been modified to 'allow' 6-node triangulations. ! However, no effort is made to actually process the midside nodes. ! Only information from the vertices is used. ! ! Modified: ! ! 07 November 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) Z(2,N), the points. ! ! Input, integer ( kind = 4 ) TRIANGLE_ORDER, the order of the triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM), ! the triangulation. ! ! Output, real ( kind = 8 ) ALPHA_MEASURE, the ALPHA quality measure. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) triangle_num integer ( kind = 4 ) triangle_order real ( kind = 8 ) a_angle integer ( kind = 4 ) a_index real ( kind = 8 ) a_x real ( kind = 8 ) a_y real ( kind = 8 ) ab_len real ( kind = 8 ) alpha real ( kind = 8 ) alpha_measure real ( kind = 8 ) arc_cosine real ( kind = 8 ) b_angle integer ( kind = 4 ) b_index real ( kind = 8 ) b_x real ( kind = 8 ) b_y real ( kind = 8 ) bc_len real ( kind = 8 ) c_angle integer ( kind = 4 ) c_index real ( kind = 8 ) c_x real ( kind = 8 ) c_y real ( kind = 8 ) ca_len real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) real ( kind = 8 ) z(2,n) alpha = huge ( alpha ) do triangle = 1, triangle_num a_index = triangle_node(1,triangle) b_index = triangle_node(2,triangle) c_index = triangle_node(3,triangle) a_x = z(1,a_index) a_y = z(2,a_index) b_x = z(1,b_index) b_y = z(2,b_index) c_x = z(1,c_index) c_y = z(2,c_index) ab_len = sqrt ( ( a_x - b_x )**2 + ( a_y - b_y )**2 ) bc_len = sqrt ( ( b_x - c_x )**2 + ( b_y - c_y )**2 ) ca_len = sqrt ( ( c_x - a_x )**2 + ( c_y - a_y )**2 ) ! ! Take care of a ridiculous special case. ! if ( ab_len == 0.0D+00 .and. & bc_len == 0.0D+00 .and. & ca_len == 0.0D+00 ) then a_angle = 2.0D+00 * pi / 3.0D+00 b_angle = 2.0D+00 * pi / 3.0D+00 c_angle = 2.0D+00 * pi / 3.0D+00 else if ( ca_len == 0.0D+00 .or. ab_len == 0.0D+00 ) then a_angle = pi else a_angle = arc_cosine ( ( ca_len**2 + ab_len**2 - bc_len**2 ) & / ( 2.0D+00 * ca_len * ab_len ) ) end if if ( ab_len == 0.0D+00 .or. bc_len == 0.0D+00 ) then b_angle = pi else b_angle = arc_cosine ( ( ab_len**2 + bc_len**2 - ca_len**2 ) & / ( 2.0D+00 * ab_len * bc_len ) ) end if if ( bc_len == 0.0D+00 .or. ca_len == 0.0D+00 ) then c_angle = pi else c_angle = arc_cosine ( ( bc_len**2 + ca_len**2 - ab_len**2 ) & / ( 2.0D+00 * bc_len * ca_len ) ) end if end if alpha = min ( alpha, a_angle ) alpha = min ( alpha, b_angle ) alpha = min ( alpha, c_angle ) end do ! ! Normalize angles from [0,60] degrees into qualities in [0,1]. ! alpha_measure = alpha * 3.0D+00 / pi return end function angle_rad_2d ( p1, p2, p3 ) !*****************************************************************************80 ! !! ANGLE_RAD_2D returns the angle 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 ! ! 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 ( p(1) == 0.0D+00 .and. p(2) == 0.0D+00 ) then angle_rad_2d = 0.0D+00 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 function arc_cosine ( c ) !*****************************************************************************80 ! !! ARC_COSINE computes the arc cosine function, with argument truncation. ! ! Discussion: ! ! If you call your system ACOS routine with an input argument that is ! even slightly outside the range [-1.0, 1.0 ], you may get an unpleasant ! surprise (I did). ! ! This routine simply truncates arguments outside the range. ! ! Modified: ! ! 02 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) C, the argument. ! ! Output, real ( kind = 8 ) ARC_COSINE, an angle whose cosine is C. ! implicit none real ( kind = 8 ) arc_cosine real ( kind = 8 ) c real ( kind = 8 ) c2 c2 = c c2 = max ( c2, -1.0D+00 ) c2 = min ( c2, +1.0D+00 ) arc_cosine = acos ( c2 ) return end function area_measure ( n, z, triangle_order, triangle_num, triangle_node ) !*****************************************************************************80 ! !! AREA_MEASURE determines the area ratio quality measure. ! ! Discusion: ! ! This measure computes the area of every triangle, and returns ! the ratio of the minimum to the maximum triangle. A value of ! 1 is "perfect", indicating that all triangles have the same area. ! A value of 0 is the worst possible result. ! ! The code has been modified to 'allow' 6-node triangulations. ! However, no effort is made to actually process the midside nodes. ! Only information from the vertices is used. ! ! Modified: ! ! 07 November 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) Z(2,N), the points. ! ! Input, integer ( kind = 4 ) TRIANGLE_ORDER, the order of the triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM), ! the triangulation. ! ! Output, real ( kind = 8 ) AREA_MEASURE, the AREA quality measure. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) triangle_num integer ( kind = 4 ) triangle_order real ( kind = 8 ) area real ( kind = 8 ) area_max real ( kind = 8 ) area_measure real ( kind = 8 ) area_min integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) z(2,n) area_max = 0.0D+00 area_min = huge ( area_min ) do triangle = 1, triangle_num x1 = z(1,triangle_node(1,triangle)) y1 = z(2,triangle_node(1,triangle)) x2 = z(1,triangle_node(2,triangle)) y2 = z(2,triangle_node(2,triangle)) x3 = z(1,triangle_node(3,triangle)) y3 = z(2,triangle_node(3,triangle)) area = 0.5D+00 * abs ( x1 * ( y2 - y3 ) & + x2 * ( y3 - y1 ) & + x3 * ( y1 - y2 ) ) area_min = min ( area_min, area ) area_max = max ( area_max, area ) end do area_measure = area_min / area_max return end subroutine bandwidth ( element_order, element_num, element_node, ml, mu, m ) !*****************************************************************************80 ! !! BANDWIDTH determines the bandwidth associated with the finite element mesh. ! ! Discussion: ! ! The quantity computed here is the "geometric" bandwidth determined ! by the finite element mesh alone. ! ! If a single finite element variable is associated with each node ! of the mesh, and if the nodes and variables are numbered in the ! same way, then the geometric bandwidth is the same as the bandwidth ! of a typical finite element matrix. ! ! The bandwidth M is defined in terms of the lower and upper bandwidths: ! ! M = ML + 1 + MU ! ! where ! ! ML = maximum distance from any diagonal entry to a nonzero ! entry in the same row, but earlier column, ! ! MU = maximum distance from any diagonal entry to a nonzero ! entry in the same row, but later column. ! ! Because the finite element node adjacency relationship is symmetric, ! we are guaranteed that ML = MU. ! ! Modified: ! ! 02 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Output, integer ( kind = 4 ) ML, MU, the lower and upper bandwidths ! of the matrix. ! ! Output, integer ( kind = 4 ) M, the bandwidth of the matrix. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) element_order integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) global_i integer ( kind = 4 ) global_j integer ( kind = 4 ) local_i integer ( kind = 4 ) local_j integer ( kind = 4 ) m integer ( kind = 4 ) ml integer ( kind = 4 ) mu ml = 0 mu = 0 do element = 1, element_num do local_i = 1, element_order global_i = element_node(local_i,element) do local_j = 1, element_order global_j = element_node(local_j,element) mu = max ( mu, global_j - global_i ) ml = max ( ml, global_i - global_j ) end do end do end do m = ml + 1 + mu 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. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! 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 get_seed ( seed ) !*****************************************************************************80 ! !! GET_SEED returns a seed for the random number generator. ! ! Discussion: ! ! The seed depends on the current time, and ought to be (slightly) ! different every millisecond. Once the seed is obtained, a random ! number generator should be called a few times to further process ! the seed. ! ! Modified: ! ! 17 November 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) SEED, a pseudorandom seed value. ! implicit none integer ( kind = 4 ) seed real ( kind = 8 ) temp character ( len = 10 ) time character ( len = 8 ) today integer ( kind = 4 ) values(8) character ( len = 5 ) zone call date_and_time ( today, time, zone, values ) temp = 0.0D+00 temp = temp + real ( values(2) - 1, kind = 8 ) / 11.0D+00 temp = temp + real ( values(3) - 1, kind = 8 ) / 30.0D+00 temp = temp + real ( values(5), kind = 8 ) / 23.0D+00 temp = temp + real ( values(6), kind = 8 ) / 59.0D+00 temp = temp + real ( values(7), kind = 8 ) / 59.0D+00 temp = temp + real ( values(8), kind = 8 ) / 999.0D+00 temp = temp / 6.0D+00 ! ! Force 0 < TEMP <= 1. ! do while ( temp <= 0.0D+00 ) temp = temp + 1.0D+00 end do do while ( 1.0D+00 < temp ) temp = temp - 1.0D+00 end do seed = int ( real ( huge ( 1 ), kind = 8 ) * temp ) ! ! Never use a seed of 0 or maximum integer. ! if ( seed == 0 ) then seed = 1 end if if ( seed == huge ( 1 ) ) then seed = seed - 1 end if return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is an integer between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) IUNIT, the free unit number. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of I4 division. ! ! Discussion: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! The MOD function computes a result with the same sign as the ! quantity being divided. Thus, suppose you had an angle A, ! and you wanted to ensure that it was between 0 and 360. ! Then mod(A,360) would do, if A was positive, but if A ! was negative, your result would be between -360 and 0. ! ! On the other hand, I4_MODP(A,360) is between 0 and 360, always. ! ! Examples: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) I, the number to be divided. ! ! Input, integer ( kind = 4 ) J, the number that divides I. ! ! Output, integer ( kind = 4 ) I4_MODP, the nonnegative remainder when I is ! divided by 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 i4_swap ( i, j ) !*****************************************************************************80 ! !! I4_SWAP switches two I4's. ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer ( kind = 4 ) I, J. On output, the values of I and ! J have been interchanged. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k k = i i = j j = k return end function i4_uniform ( a, b, seed ) !*****************************************************************************80 ! !! I4_UNIFORM returns a scaled pseudorandom I4. ! ! Discussion: ! ! An I4 is an integer ( kind = 4 ) value. ! ! The pseudorandom number will be scaled to be uniformly distributed ! between A and B. ! ! Modified: ! ! 12 November 2006 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley Interscience, page 95, 1998. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer ( kind = 4 ) A, B, the limits of the interval. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, integer ( kind = 4 ) I4_UNIFORM, a number between A and B. ! implicit none integer ( kind = 4 ) a integer ( kind = 4 ) b integer ( kind = 4 ) i4_uniform integer ( kind = 4 ) k real ( kind = 4 ) r integer ( kind = 4 ) seed integer ( kind = 4 ) value if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_UNIFORM - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if r = real ( seed, kind = 4 ) * 4.656612875E-10 ! ! Scale R to lie between A-0.5 and B+0.5. ! r = ( 1.0E+00 - r ) * ( real ( min ( a, b ), kind = 4 ) - 0.5E+00 ) & + r * ( real ( max ( a, b ), kind = 4 ) + 0.5E+00 ) ! ! Use rounding to convert R to an integer between A and B. ! value = nint ( r, kind = 4 ) value = max ( value, min ( a, b ) ) value = min ( value, max ( a, b ) ) i4_uniform = value 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 ! ! Modified: ! ! 19 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) IVAL, an integer value. ! ! Input, integer ( kind = 4 ) ILO, IHI, the desired bounds for the ! integer value. ! ! 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 ) jhi integer ( kind = 4 ) jlo integer ( kind = 4 ) wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then i4_wrap = jlo else i4_wrap = jlo + i4_modp ( ival - jlo, wide ) end if return end subroutine i4col_compare ( m, n, a, i, j, isgn ) !*****************************************************************************80 ! !! I4COL_COMPARE compares columns I and J of an I4COL. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! ISGN = -1 ! ! Modified: ! ! 30 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), an array of N columns of vectors ! of length M. ! ! Input, integer ( kind = 4 ) I, J, the columns to be compared. ! I and J must be between 1 and N. ! ! Output, integer ( kind = 4 ) ISGN, the results of the comparison: ! -1, column I < column J, ! 0, column I = column J, ! +1, column J < column I. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) i integer ( kind = 4 ) isgn integer ( kind = 4 ) j integer ( kind = 4 ) k ! ! Check. ! if ( i < 1 .or. n < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index I is out of bounds.' stop end if if ( j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index J is out of bounds.' stop end if isgn = 0 if ( i == j ) then return end if k = 1 do while ( k <= m ) if ( a(k,i) < a(k,j) ) then isgn = -1 return else if ( a(k,j) < a(k,i) ) then isgn = +1 return end if k = k + 1 end do return end subroutine i4col_sort_a ( m, n, a ) !*****************************************************************************80 ! !! I4COL_SORT_A ascending sorts an I4COL. ! ! Discussion: ! ! In lexicographic order, the statement "X < Y", applied to two real ! vectors X and Y of length M, means that there is some index I, with ! 1 <= I <= M, with the property that ! ! X(J) = Y(J) for J < I, ! and ! X(I) < Y(I). ! ! In other words, the first time they differ, X is smaller. ! ! Modified: ! ! 25 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows of A, and the length of ! a vector of data. ! ! Input, integer ( kind = 4 ) N, the number of columns of A. ! ! Input/output, integer ( kind = 4 ) A(M,N). ! On input, the array of N columns of M-vectors. ! On output, the columns of A have been sorted in ascending ! lexicographic order. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) i integer ( kind = 4 ) indx integer ( kind = 4 ) isgn integer ( kind = 4 ) j if ( m <= 0 ) then return end if if ( n <= 1 ) then return end if ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then call i4col_swap ( m, n, a, i, j ) ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4col_compare ( m, n, a, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4col_sorted_unique_count ( m, n, a, unique_num ) !*****************************************************************************80 ! !! I4COL_SORTED_UNIQUE_COUNT counts unique elements in an I4COL. ! ! Discussion: ! ! The columns of the array may be ascending or descending sorted. ! ! Modified: ! ! 17 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), a sorted array, containing ! N columns of data. ! ! Output, integer ( kind = 4 ) UNIQUE_NUM, the number of unique columns. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) j1 integer ( kind = 4 ) j2 integer ( kind = 4 ) unique_num if ( n <= 0 ) then unique_num = 0 return end if unique_num = 1 j1 = 1 do j2 = 2, n if ( any ( a(1:m,j1) /= a(1:m,j2) ) ) then unique_num = unique_num + 1 j1 = j2 end if end do return end subroutine i4col_swap ( m, n, a, i, j ) !*****************************************************************************80 ! !! I4COL_SWAP swaps columns I and J of an I4COL. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! A = ( ! 1 4 3 2 ! 5 8 7 6 ! 9 12 11 10 ) ! ! Modified: ! ! 04 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns ! in the array. ! ! Input/output, integer ( kind = 4 ) A(M,N), an array of N columns ! of length M. ! ! Input, integer ( kind = 4 ) I, J, the columns to be swapped. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) col(m) integer ( kind = 4 ) i integer ( kind = 4 ) j if ( i < 1 .or. n < i .or. j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_SWAP - Fatal error!' write ( *, '(a)' ) ' I or J is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i write ( *, '(a,i8)' ) ' J = ', j write ( *, '(a,i8)' ) ' N = ', n stop end if if ( i == j ) then return end if col(1:m) = a(1:m,i) a(1:m,i) = a(1:m,j) a(1:m,j) = col(1:m) return end subroutine i4i4_sort_a ( i1, i2, j1, j2 ) !*****************************************************************************80 ! !! I4I4_SORT_A ascending sorts a pair of integers. ! ! Discussion: ! ! An I4I4 is a pair of integers, regarded as a single data item. ! ! The program allows the reasonable call: ! ! call i4i4_sort_a ( i1, i2, i1, i2 ) ! ! and this will return the reasonable result. ! ! Modified: ! ! 11 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) I1, I2, the values to sort. ! ! Output, integer ( kind = 4 ) J1, J2, the sorted values. ! implicit none integer ( kind = 4 ) i1 integer ( kind = 4 ) i2 integer ( kind = 4 ) j1 integer ( kind = 4 ) j2 integer ( kind = 4 ) k1 integer ( kind = 4 ) k2 ! ! Copy arguments, so that the user can make "reasonable" calls like: ! ! call i4i4_sort_a ( i1, i2, i1, i2 ) ! k1 = i1 k2 = i2 j1 = min ( k1, k2 ) j2 = max ( k1, k2 ) return end subroutine i4mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed. ! ! Modified: ! ! 28 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ), parameter :: i4_1 = 1 character ( len = * ) title call i4mat_transpose_print_some ( m, n, a, i4_1, i4_1, m, n, title ) return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT. ! ! Modified: ! ! 09 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ), parameter :: incx = 10 integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) character ( len = 7 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2 integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine i4vec_heap_d ( n, a ) !*****************************************************************************80 ! !! I4VEC_HEAP_D reorders an I4VEC into a descending heap. ! ! Discussion: ! ! A descending heap is an array A with the property that, for every index J, ! A(2*J) <= A(J) and A(2*J+1) <= A(J), (as long as the indices ! 2*J and 2*J+1 are legal). ! ! Diagram: ! ! A(1) ! / \ ! A(2) A(3) ! / \ / \ ! A(4) A(5) A(6) A(7) ! / \ / \ ! A(8) A(9) A(10) A(11) ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the size of the input array. ! ! Input/output, integer ( kind = 4 ) A(N). ! On input, an unsorted array. ! On output, the array has been reordered into a heap. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) i integer ( kind = 4 ) ifree integer ( kind = 4 ) key integer ( kind = 4 ) m ! ! Only nodes N/2 down to 1 can be "parent" nodes. ! do i = n/2, 1, -1 ! ! Copy the value out of the parent node. ! Position IFREE is now "open". ! key = a(i) ifree = i do ! ! Positions 2*IFREE and 2*IFREE + 1 are the descendants of position ! IFREE. (One or both may not exist because they exceed N.) ! m = 2 * ifree ! ! Does the first position exist? ! if ( n < m ) then exit end if ! ! Does the second position exist? ! if ( m + 1 <= n ) then ! ! If both positions exist, take the larger of the two values, ! and update M if necessary. ! if ( a(m) < a(m+1) ) then m = m + 1 end if end if ! ! If the large descendant is larger than KEY, move it up, ! and update IFREE, the location of the free position, and ! consider the descendants of THIS position. ! if ( a(m) <= key ) then exit end if a(ifree) = a(m) ifree = m end do ! ! Once there is no more shifting to do, KEY moves into the free spot IFREE. ! a(ifree) = key end do return end subroutine i4vec_indicator ( n, a ) !*****************************************************************************80 ! !! I4VEC_INDICATOR sets an I4VEC to the indicator vector. ! ! 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 i4vec_print ( n, a, title ) !*****************************************************************************80 ! !! I4VEC_PRINT prints an I4VEC. ! ! Modified: ! ! 28 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of components of the vector. ! ! Input, integer ( kind = 4 ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) big integer ( kind = 4 ) i character ( len = * ) title if ( title /= ' ' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if big = maxval ( abs ( a(1:n) ) ) write ( *, '(a)' ) ' ' if ( big < 1000 ) then do i = 1, n write ( *, '(i8,1x,i4)' ) i, a(i) end do else if ( big < 1000000 ) then do i = 1, n write ( *, '(i8,1x,i7)' ) i, a(i) end do else do i = 1, n write ( *, '(i8,i11)' ) i, a(i) end do end if return end subroutine i4vec_sort_heap_a ( n, a ) !*****************************************************************************80 ! !! I4VEC_SORT_HEAP_A ascending sorts an I4VEC using heap sort. ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the array. ! ! Input/output, integer ( kind = 4 ) A(N). ! On input, the array to be sorted; ! On output, the array has been sorted. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) n1 if ( n <= 1 ) then return end if ! ! 1: Put A into descending heap form. ! call i4vec_heap_d ( n, a ) ! ! 2: Sort A. ! ! The largest object in the heap is in A(1). ! Move it to position A(N). ! call i4_swap ( a(1), a(n) ) ! ! Consider the diminished heap of size N1. ! do n1 = n-1, 2, -1 ! ! Restore the heap structure of A(1) through A(N1). ! call i4vec_heap_d ( n1, a ) ! ! Take the largest object from A(1) and move it to A(N1). ! call i4_swap ( a(1), a(n1) ) end do return end subroutine i4vec_sorted_unique ( n, a, unique_num ) !*****************************************************************************80 ! !! I4VEC_SORTED_UNIQUE gets the unique elements in a sorted I4VEC. ! ! Modified: ! ! 09 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of elements in A. ! ! Input/output, integer ( kind = 4 ) A(N). On input, the sorted ! integer array. On output, the unique elements in A. ! ! Output, integer ( kind = 4 ) UNIQUE_NUM, the number of unique elements ! in A. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) itest integer ( kind = 4 ) unique_num unique_num = 0 if ( n <= 0 ) then return end if unique_num = 1 do itest = 2, n if ( a(itest) /= a(unique_num) ) then unique_num = unique_num + 1 a(unique_num) = a(itest) end if end do return end subroutine i4vec2_compare ( n, a1, a2, i, j, isgn ) !*****************************************************************************80 ! !! I4VEC2_COMPARE compares pairs of integers stored in two vectors. ! ! Modified: ! ! 22 October 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of data items. ! ! Input, integer ( kind = 4 ) A1(N), A2(N), contain the two components ! of each item. ! ! Input, integer ( kind = 4 ) I, J, the items to be compared. ! ! Output, integer ( kind = 4 ) ISGN, the results of the comparison: ! -1, item I < item J, ! 0, item I = item J, ! +1, item J < item I. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a1(n) integer ( kind = 4 ) a2(n) integer ( kind = 4 ) i integer ( kind = 4 ) isgn integer ( kind = 4 ) j isgn = 0 if ( a1(i) < a1(j) ) then isgn = -1 else if ( a1(i) == a1(j) ) then if ( a2(i) < a2(j) ) then isgn = -1 else if ( a2(i) < a2(j) ) then isgn = 0 else if ( a2(j) < a2(i) ) then isgn = +1 end if else if ( a1(j) < a1(i) ) then isgn = +1 end if return end subroutine i4vec2_sort_a ( n, a1, a2 ) !*****************************************************************************80 ! !! I4VEC2_SORT_A ascending sorts a vector of pairs of integers. ! ! Discussion: ! ! Each item to be sorted is a pair of integers (I,J), with the I ! and J values stored in separate vectors A1 and A2. ! ! Modified: ! ! 25 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of items of data. ! ! Input/output, integer ( kind = 4 ) A1(N), A2(N), the data to be sorted. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a1(n) integer ( kind = 4 ) a2(n) integer ( kind = 4 ) i integer ( kind = 4 ) indx integer ( kind = 4 ) isgn integer ( kind = 4 ) j integer ( kind = 4 ) temp if ( n <= 1 ) then return end if ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then temp = a1(i) a1(i) = a1(j) a1(j) = temp temp = a2(i) a2(i) = a2(j) a2(j) = temp ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4vec2_compare ( n, a1, a2, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4vec2_sorted_unique ( n, a1, a2, unique_num ) !*****************************************************************************80 ! !! I4VEC2_SORTED_UNIQUE gets the unique elements in a sorted I4VEC2. ! ! Discussion: ! ! Item I is stored as the pair A1(I), A2(I). ! ! The items must have been sorted, or at least it must be the ! case that equal items are stored in adjacent vector locations. ! ! If the items were not sorted, then this routine will only ! replace a string of equal values by a single representative. ! ! Modified: ! ! 09 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of items. ! ! Input/output, integer ( kind = 4 ) A1(N), A2(N). ! On input, the array of N items. ! On output, an array of unique items. ! ! Output, integer ( kind = 4 ) UNIQUE_NUM, the number of unique items. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a1(n) integer ( kind = 4 ) a2(n) integer ( kind = 4 ) itest integer ( kind = 4 ) unique_num unique_num = 0 if ( n <= 0 ) then return end if unique_num = 1 do itest = 2, n if ( a1(itest) /= a1(unique_num) .or. a2(itest) /= a2(unique_num) ) then unique_num = unique_num + 1 a1(unique_num) = a1(itest) a2(unique_num) = a2(itest) end if end do 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). ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! 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 lvec_print ( n, a, title ) !*****************************************************************************80 ! !! LVEC_PRINT prints an LVEC. ! ! Modified: ! ! 26 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of components of the vector. ! ! Input, logical A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer ( kind = 4 ) n logical a(n) integer ( kind = 4 ) i character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,l1)' ) i, a(i) end do return end subroutine node_merge ( dim_num, node_num, node_xy, tolerance, node_rep ) !*****************************************************************************80 ! !! NODE_MERGE detects nodes that should be merged. ! ! Discussion: ! ! Two nodes "should" be merged if they are within TOLERANCE distance ! of each other. ! ! With a tolerance of 0, only exactly equal nodes are counted. ! ! With a positive tolerance, a pair of nodes inside a circle of ! radius TOLERANCE result in a count of 1 duplicate. ! ! However, what do we do if nodes A, B and C are arranged in a line,! ! with A and B just within TOLERANCE of each other, and B and C just ! within tolerance of each other? What we do here is make a choice ! that can be defended consistently. A and B define an equivalence ! class because they are closer than TOLERANCE. C is then added to ! this equivalence class, because it is within TOLERANCE of at least ! on thing in that equivalence class. ! ! Thus, if 100 nodes are separated pairwise by slightly less ! than TOLERANCE, a total of 99 duplicates will be counted. ! ! The program starts out by giving each node its own label. ! If it finds that two nodes should be merged, then the index of ! one node is used as the label for both. This process continues ! until all nodes have been considered. The number of unique nodes ! is the number of unique values in the output quantity NODE_REP. ! ! Modified: ! ! 14 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension. ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(DIM_NUM,NODE_NUM), the nodes. ! ! Input, real ( kind = 8 ) TOLERANCE, the maximum distance between ! two nodes regarded as duplicate. ! ! Output, integer ( kind = 4 ) NODE_REP(NODE_NUM), the "representative" of ! each node. NODE_REP(NODE) is the index of a node which is within ! TOLERANCE of node NODE, or for which a chain of nodes can be found, all ! having the same representative, and all of which are pairwise closer ! than TOLERANCE. ! implicit none integer ( kind = 4 ) dim_num integer ( kind = 4 ) node_num real ( kind = 8 ) dist integer ( kind = 4 ) node_rep(node_num) real ( kind = 8 ) node_xy(dim_num,node_num) integer ( kind = 4 ) node1 integer ( kind = 4 ) node2 integer ( kind = 4 ) rep real ( kind = 8 ) rep_dist(node_num) real ( kind = 8 ) tolerance do node1 = 1, node_num node_rep(node1) = node1 end do do node1 = 1, node_num rep_dist(1:node_num) = huge ( 1.0D+00 ) do node2 = 1, node_num dist = sqrt ( sum ( & ( node_xy(1:dim_num,node1) - node_xy(1:dim_num,node2) )**2 ) ) rep = node_rep(node2) if ( dist < rep_dist(rep) ) then rep_dist(rep) = dist end if end do do node2 = 1, node_num rep = node_rep(node2) if ( rep_dist(rep) <= tolerance ) then node_rep(node2) = node1 end if end do end do return end subroutine ns_adj_col_set ( node_num, triangle_num, variable_num, & triangle_node, triangle_neighbor, node_u_variable, node_v_variable, & node_p_variable, adj_num, adj_col ) !*****************************************************************************80 ! !! NS_ADJ_COL_SET sets the COL array in a Navier Stokes triangulation. ! ! Discussion: ! ! This routine also returns the value of ADJ_NUM, the number of ! Navier-Stokes variable adjacencies, which should be identical to the ! value that would have been computed by calling NS_ADJ_COUNT. ! ! This routine is called to set up the ADJ_COL array, which indicates ! the number of entries needed to store each column in the sparse ! compressed column storage to be used for the adjancency matrix. ! ! The triangulation is assumed to involve 6-node triangles. ! ! Variables for the horizontal and vertical velocities are associated ! with every node. Variables for the pressure are associated only with ! the vertex nodes. ! ! We are interested in determining the number of nonzero entries in the ! stiffness matrix of the Stokes equations, or the jacobian matrix of ! the Navier Stokes equations. To this end, we will say, somewhat ! too broadly, that two variables are "adjacent" if their associated ! nodes both occur in some common element. This adjacency of variables ! I and J is taken to be equivalent to the possible nonzeroness of ! matrix entries A(I,J) and A(J,I). ! ! A sparse compressed column format is used to store the counts for ! the nonzeroes. In other words, while the value ADJ_NUM reports the ! number of adjacencies, the vector ADJ_COL is sufficient to allow us ! to properly set up a sparse compressed matrix for the actual storage ! of the sparse matrix, if we desire to proceed. ! ! Local Node Numbering: ! ! 3 ! s |\ ! i | \ ! d | \ ! e 6 5 side 2 ! | \ ! 3 | \ ! | \ ! 1---4---2 ! ! side 1 ! ! Variable Diagram: ! ! UVP ! |\ ! | \ ! | \ ! UV UV ! | \ ! | \ ! | \ ! UVP--UV--UVP ! ! Modified: ! ! 26 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) VARIABLE_NUM, the number of variables. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(6,TRIANGLE_NUM), lists the ! nodes that make up each triangle. The first three nodes are the vertices, ! in counterclockwise order. The fourth value is the midside ! node between nodes 1 and 2; the fifth and sixth values are ! the other midside nodes in the logical order. ! ! Input, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), for each ! side of a triangle, lists the neighboring triangle, or -1 if there is ! no neighbor. ! ! Input, integer ( kind = 4 ) NODE_U_VARIABLE(NODE_NUM), ! NODE_V_VARIABLE(NODE_NUM), NODE_P_VARIABLE(NODE_NUM), the index of the ! horizontal velocity, vertical velocity and pressure variables associated ! with a node, or -1 if no such variable is associated with the node. ! ! Output, integer ( kind = 4 ) ADJ_NUM, the number of ! Navier Stokes variable adjacencies. ! ! Output, integer ( kind = 4 ) ADJ_COL(VARIABLE_NUM+1). Information about ! variable J is stored in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 6 integer ( kind = 4 ) variable_num integer ( kind = 4 ) adj_num integer ( kind = 4 ) adj_col(variable_num+1) integer ( kind = 4 ) n1 integer ( kind = 4 ) n2 integer ( kind = 4 ) n3 integer ( kind = 4 ) n4 integer ( kind = 4 ) n5 integer ( kind = 4 ) n6 integer ( kind = 4 ) node integer ( kind = 4 ) node_p_variable(node_num) integer ( kind = 4 ) node_u_variable(node_num) integer ( kind = 4 ) node_v_variable(node_num) integer ( kind = 4 ) p1 integer ( kind = 4 ) p2 integer ( kind = 4 ) p3 integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle2 integer ( kind = 4 ) triangle_neighbor(3,triangle_num) integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) integer ( kind = 4 ) u1 integer ( kind = 4 ) u2 integer ( kind = 4 ) u3 integer ( kind = 4 ) u4 integer ( kind = 4 ) u5 integer ( kind = 4 ) u6 integer ( kind = 4 ) v1 integer ( kind = 4 ) v2 integer ( kind = 4 ) v3 integer ( kind = 4 ) v4 integer ( kind = 4 ) v5 integer ( kind = 4 ) v6 integer ( kind = 4 ) variable adj_num = 0 ! ! Set every variable to be adjacent to itself. ! adj_col(1:variable_num) = 1 ! ! Set every variable to be adjacent to the other variables associated with ! that node. ! ! U <=> V ! U <=> P (if there is a P variable) ! V <=> P (if there is a P variable) ! do node = 1, node_num u1 = node_u_variable(node) v1 = node_v_variable(node) p1 = node_p_variable(node) adj_col(u1) = adj_col(u1) + 1 adj_col(v1) = adj_col(v1) + 1 if ( 0 < p1 ) then adj_col(u1) = adj_col(u1) + 1 adj_col(v1) = adj_col(v1) + 1 adj_col(p1) = adj_col(p1) + 2 end if end do ! ! Examine each triangle. ! do triangle = 1, triangle_num n1 = triangle_node(1,triangle) n2 = triangle_node(2,triangle) n3 = triangle_node(3,triangle) n4 = triangle_node(4,triangle) n5 = triangle_node(5,triangle) n6 = triangle_node(6,triangle) u1 = node_u_variable(n1) v1 = node_v_variable(n1) p1 = node_p_variable(n1) u2 = node_u_variable(n2) v2 = node_v_variable(n2) p2 = node_p_variable(n2) u3 = node_u_variable(n3) v3 = node_v_variable(n3) p3 = node_p_variable(n3) u4 = node_u_variable(n4) v4 = node_v_variable(n4) u5 = node_u_variable(n5) v5 = node_v_variable(n5) u6 = node_u_variable(n6) v6 = node_v_variable(n6) ! ! For sure, we add the new adjacencies: ! ! U5 V5 <=> U1 V1 P1 ! U6 V6 <=> U2 V2 P2 ! U4 V4 <=> U3 V3 P3 ! U5 V5 <=> U4 V4 ! U6 V6 <=> U4 V4 ! U6 V6 <=> U5 V5 ! adj_col(u1) = adj_col(u1) + 2 adj_col(v1) = adj_col(v1) + 2 adj_col(p1) = adj_col(p1) + 2 adj_col(u2) = adj_col(u2) + 2 adj_col(v2) = adj_col(v2) + 2 adj_col(p2) = adj_col(p2) + 2 adj_col(u3) = adj_col(u3) + 2 adj_col(v3) = adj_col(v3) + 2 adj_col(p3) = adj_col(p3) + 2 adj_col(u4) = adj_col(u4) + 7 adj_col(v4) = adj_col(v4) + 7 adj_col(u5) = adj_col(u5) + 7 adj_col(v5) = adj_col(v5) + 7 adj_col(u6) = adj_col(u6) + 7 adj_col(v6) = adj_col(v6) + 7 ! ! Add edges (1,2), (1,4), (2,4) if this is the first occurrence, ! that is, if the edge (1,4,2) is on a boundary (TRIANGLE2 <= 0) ! or if this triangle is the first of the pair in which the edge ! occurs (TRIANGLE < TRIANGLE2). ! ! Maybe add ! ! U1 V1 P1 <=> U2 V2 P2 ! U1 V1 P1 <=> U4 V4 ! U2 V2 P2 <=> U4 V4 ! triangle2 = triangle_neighbor(1,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_col(u1) = adj_col(u1) + 5 adj_col(v1) = adj_col(v1) + 5 adj_col(p1) = adj_col(p1) + 5 adj_col(u2) = adj_col(u2) + 5 adj_col(v2) = adj_col(v2) + 5 adj_col(p2) = adj_col(p2) + 5 adj_col(u4) = adj_col(u4) + 6 adj_col(v4) = adj_col(v4) + 6 end if ! ! Maybe add ! ! U2 V2 P2 <=> U3 V3 P3 ! U2 V2 P2 <=> U5 V5 ! U3 V3 P3 <=> U5 V5 ! triangle2 = triangle_neighbor(2,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_col(u2) = adj_col(u2) + 5 adj_col(v2) = adj_col(v2) + 5 adj_col(p2) = adj_col(p2) + 5 adj_col(u3) = adj_col(u3) + 5 adj_col(v3) = adj_col(v3) + 5 adj_col(p3) = adj_col(p3) + 5 adj_col(u5) = adj_col(u5) + 6 adj_col(v5) = adj_col(v5) + 6 end if ! ! Maybe add ! ! U1 V1 P1 <=> U3 V3 P3 ! U1 V1 P1 <=> U6 V6 ! U3 V3 P3 <=> U6 V6 ! triangle2 = triangle_neighbor(3,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_col(u1) = adj_col(u1) + 5 adj_col(v1) = adj_col(v1) + 5 adj_col(p1) = adj_col(p1) + 5 adj_col(u3) = adj_col(u3) + 5 adj_col(v3) = adj_col(v3) + 5 adj_col(p3) = adj_col(p3) + 5 adj_col(u6) = adj_col(u6) + 6 adj_col(v6) = adj_col(v6) + 6 end if end do ! ! We used ADJ_COL to count the number of entries in each column. ! Convert it to pointers into the ADJ array. ! adj_col(2:variable_num+1) = adj_col(1:variable_num) adj_col(1) = 1 do variable = 2, variable_num + 1 adj_col(variable) = adj_col(variable-1) + adj_col(variable) end do adj_num = adj_col(variable_num+1) - 1 return end subroutine ns_adj_count ( node_num, triangle_num, variable_num, triangle_node, & triangle_neighbor, node_u_variable, node_v_variable, node_p_variable, & adj_num ) !*****************************************************************************80 ! !! NS_ADJ_COUNT counts adjacencies in a Navier Stokes triangulation. ! ! Discussion: ! ! This routine is called to count the adjacencies, so that the ! appropriate amount of memory can be set aside for storage when ! the adjacency structure is created. ! ! The value of ADJ_NUM computed and returned by this routine should ! be identical to the value computed by NS_ADJ_COL_SET. ! ! The triangulation is assumed to involve 6-node triangles. ! ! Variables for the horizontal and vertical velocities are associated ! with every node. Variables for the pressure are associated only with ! the vertex nodes. ! ! We are interested in determining the number of nonzero entries in the ! stiffness matrix of the Stokes equations, or the jacobian matrix of ! the Navier Stokes equations. To this end, we will say, somewhat ! too broadly, that two variables are "adjacent" if their associated ! nodes both occur in some common element. This adjacency of variables ! I and J is taken to be equivalent to the possible nonzeroness of ! matrix entries A(I,J) and A(J,I). ! ! A sparse compressed column format is used to store the counts for ! the nonzeroes. In other words, while the value ADJ_NUM reports the ! number of adjacencies, the vector ADJ_COL is sufficient to allow us ! to properly set up a sparse compressed matrix for the actual storage ! of the sparse matrix, if we desire to proceed. ! ! Local Node Numbering: ! ! 3 ! s |\ ! i | \ ! d | \ ! e 6 5 side 2 ! | \ ! 3 | \ ! | \ ! 1---4---2 ! ! side 1 ! ! Variable Diagram: ! ! UVP ! |\ ! | \ ! | \ ! UV UV ! | \ ! | \ ! | \ ! UVP--UV--UVP ! ! Modified: ! ! 19 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) VARIABLE_NUM, the number of variables. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(6,TRIANGLE_NUM), lists the ! nodes that make up each triangle. The first three nodes are the vertices, ! in counterclockwise order. The fourth value is the midside ! node between nodes 1 and 2; the fifth and sixth values are ! the other midside nodes in the logical order. ! ! Input, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), for each ! side of a triangle, lists the neighboring triangle, or -1 if there is ! no neighbor. ! ! Input, integer ( kind = 4 ) NODE_U_VARIABLE(NODE_NUM), ! NODE_V_VARIABLE(NODE_NUM), NODE_P_VARIABLE(NODE_NUM), the index of the ! horizontal velocity, vertical velocity and pressure variables associated ! with a node, or -1 if no such variable is associated with the node. ! ! Output, integer ( kind = 4 ) ADJ_NUM, the value of ADJ_NUM, the number of ! Navier Stokes variable adjacencies. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 6 integer ( kind = 4 ) variable_num integer ( kind = 4 ) adj_num integer ( kind = 4 ) node integer ( kind = 4 ) node_p_variable(node_num) integer ( kind = 4 ) node_u_variable(node_num) integer ( kind = 4 ) node_v_variable(node_num) integer ( kind = 4 ) p1 integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle2 integer ( kind = 4 ) triangle_neighbor(3,triangle_num) integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) integer ( kind = 4 ) variable adj_num = 0 ! ! Set every variable to be adjacent to itself. ! adj_num = variable_num ! ! Set every variable to be adjacent to the other variables associated with ! that node. ! ! U <=> V ! U <=> P (if there is a P variable) ! V <=> P (if there is a P variable) ! do node = 1, node_num adj_num = adj_num + 2 p1 = node_p_variable(node) if ( 0 < p1 ) then adj_num = adj_num + 4 end if end do ! ! Examine each triangle. ! do triangle = 1, triangle_num ! ! For sure, we add the new adjacencies: ! ! U5 V5 <=> U1 V1 P1 ! U6 V6 <=> U2 V2 P2 ! U4 V4 <=> U3 V3 P3 ! U5 V5 <=> U4 V4 ! U6 V6 <=> U4 V4 ! U6 V6 <=> U5 V5 ! adj_num = adj_num + 60 ! ! Add edges (1,2), (1,4), (2,4) if this is the first occurrence, ! that is, if the edge (1,4,2) is on a boundary (TRIANGLE2 <= 0) ! or if this triangle is the first of the pair in which the edge ! occurs (TRIANGLE < TRIANGLE2). ! ! Maybe add ! ! U1 V1 P1 <=> U2 V2 P2 ! U1 V1 P1 <=> U4 V4 ! U2 V2 P2 <=> U4 V4 ! triangle2 = triangle_neighbor(1,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_num = adj_num + 42 end if ! ! Maybe add ! ! U2 V2 P2 <=> U3 V3 P3 ! U2 V2 P2 <=> U5 V5 ! U3 V3 P3 <=> U5 V5 ! triangle2 = triangle_neighbor(2,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_num = adj_num + 42 end if ! ! Maybe add ! ! U1 V1 P1 <=> U3 V3 P3 ! U1 V1 P1 <=> U6 V6 ! U3 V3 P3 <=> U6 V6 ! triangle2 = triangle_neighbor(3,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_num = adj_num + 42 end if end do return end subroutine ns_adj_insert ( v1, v2, variable_num, adj_num, adj_col_free, & adj_row ) !*****************************************************************************80 ! !! NS_ADJ_INSERT inserts an adjacency into a compressed column adjacency matrix. ! ! Modified: ! ! 19 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) V1, V2, the indices of two items which are ! adjacent. ! ! Input, integer ( kind = 4 ) VARIABLE_NUM, the number of items. ! ! Input, integer ( kind = 4 ) ADJ_NUM, the number of entries available ! in ADJ_ROW. ! ! Input/output, integer ( kind = 4 ) ADJ_COL_FREE(VARIABLE_NUM), the next ! free location in which an entry for a given column can be stored. On ! output, two pointers have been updated. ! ! Input/output, integer ( kind = 4 ) ADJ_ROW(ADJ_NUM), the row indices of ! the Navier Stokes variable adjacency matrix. On output, two new entries ! have been added. ! integer ( kind = 4 ) adj_num integer ( kind = 4 ) variable_num integer ( kind = 4 ) adj_col_free(variable_num) integer ( kind = 4 ) adj_row(adj_num) integer ( kind = 4 ) v1 integer ( kind = 4 ) v2 adj_row(adj_col_free(v1)) = v2 adj_col_free(v1) = adj_col_free(v1) + 1 if ( v1 == v2 ) then return end if adj_row(adj_col_free(v2)) = v1 adj_col_free(v2) = adj_col_free(v2) + 1 return end subroutine ns_adj_row_set ( node_num, triangle_num, variable_num, & triangle_node, triangle_neighbor, node_u_variable, node_v_variable, & node_p_variable, adj_num, adj_col, adj_row ) !*****************************************************************************80 ! !! NS_ADJ_ROW_SET sets the Navier Stokes sparse compressed column row indices. ! ! Discussion: ! ! After NS_ADJ_COUNT has been called to count ADJ_NUM, the number of ! variable adjacencies and to set up ADJ_COL, the compressed column pointer, ! this routine can be called to assign values to ADJ_ROW, the row ! indices of the sparse compressed column adjacency matrix. ! ! Modified: ! ! 19 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) VARIABLE_NUM, the number of variables. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(6,TRIANGLE_NUM), lists the ! nodes that make up each triangle. The first three nodes are the vertices, ! in counterclockwise order. The fourth value is the midside ! node between nodes 1 and 2; the fifth and sixth values are ! the other midside nodes in the logical order. ! ! Input, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), for each ! side of a triangle, lists the neighboring triangle, or -1 if there is ! no neighbor. ! ! Input, integer ( kind = 4 ) NODE_U_VARIABLE(NODE_NUM), ! NODE_V_VARIABLE(NODE_NUM), NODE_P_VARIABLE(NODE_NUM), the index of the ! horizontal velocity, vertical velocity and pressure variables associated ! with a node, or -1 if no such variable is associated with the node. ! ! Input, integer ( kind = 4 ) ADJ_NUM, the number of Navier Stokes variable ! adjacencies. ! ! Input, integer ( kind = 4 ) ADJ_COL(VARIABLE_NUM+1). Information about ! variable J is stored in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ. ! ! Output, integer ( kind = 4 ) ADJ_ROW(ADJ_NUM), the row indices of the ! Navier Stokes variable adjacency matrix. ! ! Local Parameters: ! ! Local, integer ADJ_COL_FREE(VARIABLE_NUM), for each column, ! the location in ADJ_ROW which can store the next row index. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 6 integer ( kind = 4 ) variable_num integer ( kind = 4 ) adj_num integer ( kind = 4 ) adj_col(variable_num+1) integer ( kind = 4 ) adj_col_free(variable_num) integer ( kind = 4 ) adj_row(adj_num) integer ( kind = 4 ) k1 integer ( kind = 4 ) k2 integer ( kind = 4 ) n1 integer ( kind = 4 ) n2 integer ( kind = 4 ) n3 integer ( kind = 4 ) n4 integer ( kind = 4 ) n5 integer ( kind = 4 ) n6 integer ( kind = 4 ) node integer ( kind = 4 ) node_p_variable(node_num) integer ( kind = 4 ) node_u_variable(node_num) integer ( kind = 4 ) node_v_variable(node_num) integer ( kind = 4 ) number integer ( kind = 4 ) p1 integer ( kind = 4 ) p2 integer ( kind = 4 ) p3 integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle2 integer ( kind = 4 ) triangle_neighbor(3,triangle_num) integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) integer ( kind = 4 ) u1 integer ( kind = 4 ) u2 integer ( kind = 4 ) u3 integer ( kind = 4 ) u4 integer ( kind = 4 ) u5 integer ( kind = 4 ) u6 integer ( kind = 4 ) v integer ( kind = 4 ) v1 integer ( kind = 4 ) v2 integer ( kind = 4 ) v3 integer ( kind = 4 ) v4 integer ( kind = 4 ) v5 integer ( kind = 4 ) v6 adj_row(1:adj_num) = -1 adj_col_free(1:variable_num) = adj_col(1:variable_num) ! ! Set every variable to be adjacent to itself. ! do v = 1, variable_num call ns_adj_insert ( v, v, variable_num, adj_num, adj_col_free, adj_row ) end do ! ! Set every variable to be adjacent to the other variables associated with ! that node. ! ! U <=> V ! U <=> P (if there is a P variable) ! V <=> P (if there is a P variable) ! do node = 1, node_num u1 = node_u_variable(node) v1 = node_v_variable(node) p1 = node_p_variable(node) call ns_adj_insert ( u1, v1, variable_num, adj_num, adj_col_free, adj_row ) if ( 0 < p1 ) then call ns_adj_insert ( u1, p1, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, p1, variable_num, adj_num, adj_col_free, & adj_row ) end if end do ! ! Examine each triangle. ! do triangle = 1, triangle_num n1 = triangle_node(1,triangle) n2 = triangle_node(2,triangle) n3 = triangle_node(3,triangle) n4 = triangle_node(4,triangle) n5 = triangle_node(5,triangle) n6 = triangle_node(6,triangle) u1 = node_u_variable(n1) v1 = node_v_variable(n1) p1 = node_p_variable(n1) u2 = node_u_variable(n2) v2 = node_v_variable(n2) p2 = node_p_variable(n2) u3 = node_u_variable(n3) v3 = node_v_variable(n3) p3 = node_p_variable(n3) u4 = node_u_variable(n4) v4 = node_v_variable(n4) u5 = node_u_variable(n5) v5 = node_v_variable(n5) u6 = node_u_variable(n6) v6 = node_v_variable(n6) ! ! For sure, we add the new adjacencies: ! ! U5 V5 <=> U1 V1 P1 ! U6 V6 <=> U2 V2 P2 ! U4 V4 <=> U3 V3 P3 ! U5 V5 <=> U4 V4 ! U6 V6 <=> U4 V4 ! U6 V6 <=> U5 V5 ! call ns_adj_insert ( u5, u1, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u5, v1, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u5, p1, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v5, u1, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v5, v1, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v5, p1, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, u2, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, v2, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, p2, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, u2, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, v2, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, p2, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u4, u3, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u4, v3, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u4, p3, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v4, u3, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v4, v3, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v4, p3, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u5, u4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u5, v4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v5, u4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v5, v4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, u4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, v4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, u4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, v4, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, u5, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( u6, v5, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, u5, variable_num, adj_num, adj_col_free, adj_row ) call ns_adj_insert ( v6, v5, variable_num, adj_num, adj_col_free, adj_row ) ! ! Add edges (1,2), (1,4), (2,4) if this is the first occurrence, ! that is, if the edge (1,4,2) is on a boundary (TRIANGLE2 <= 0) ! or if this triangle is the first of the pair in which the edge ! occurs (TRIANGLE < TRIANGLE2). ! ! Maybe add ! ! U1 V1 P1 <=> U2 V2 P2 ! U1 V1 P1 <=> U4 V4 ! U2 V2 P2 <=> U4 V4 ! triangle2 = triangle_neighbor(1,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then call ns_adj_insert ( u1, u2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, v2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, p2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, u2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, v2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, p2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, u2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, v2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, p2, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, u4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, v4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, u4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, v4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, u4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, v4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u2, u4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u2, v4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, u4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, v4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, u4, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, v4, variable_num, adj_num, adj_col_free, & adj_row ) end if ! ! Maybe add ! ! U2 V2 P2 <=> U3 V3 P3 ! U2 V2 P2 <=> U5 V5 ! U3 V3 P3 <=> U5 V5 ! triangle2 = triangle_neighbor(2,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then call ns_adj_insert ( u2, u3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u2, v3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u2, p3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, u3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, v3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, p3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, u3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, v3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, p3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u2, u5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u2, v5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, u5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v2, v5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, u5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p2, v5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u3, u5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u3, v5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v3, u5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v3, v5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p3, u5, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p3, v5, variable_num, adj_num, adj_col_free, & adj_row ) end if ! ! Maybe add ! ! U1 V1 P1 <=> U3 V3 P3 ! U1 V1 P1 <=> U6 V6 ! U3 V3 P3 <=> U6 V6 ! triangle2 = triangle_neighbor(3,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then call ns_adj_insert ( u1, u3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, v3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, p3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, u3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, v3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, p3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, u3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, v3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, p3, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, u6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u1, v6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, u6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v1, v6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, u6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p1, v6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u3, u6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( u3, v6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v3, u6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( v3, v6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p3, u6, variable_num, adj_num, adj_col_free, & adj_row ) call ns_adj_insert ( p3, v6, variable_num, adj_num, adj_col_free, & adj_row ) end if end do ! ! Ascending sort the entries for each variable. ! do v = 1, variable_num k1 = adj_col(v) k2 = adj_col(v+1)-1 number = k2 + 1 - k1 call i4vec_sort_heap_a ( number, adj_row(k1:k2) ) end do return end subroutine perm_check ( n, p, ierror ) !*****************************************************************************80 ! !! PERM_CHECK checks that a vector represents a permutation. ! ! Discussion: ! ! The routine verifies that each of the integers from 1 ! to N occurs among the N entries of the permutation. ! ! Modified: ! ! 01 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries. ! ! Input, integer ( kind = 4 ) P(N), the array to check. ! ! Output, integer ( kind = 4 ) IERROR, error flag. ! 0, the array represents a permutation. ! nonzero, the array does not represent a permutation. The smallest ! missing value is equal to IERROR. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) ierror integer ( kind = 4 ) ifind integer ( kind = 4 ) iseek integer ( kind = 4 ) p(n) ierror = 0 do iseek = 1, n ierror = iseek do ifind = 1, n if ( p(ifind) == iseek ) then ierror = 0 exit end if end do if ( ierror /= 0 ) then return end if end do return end subroutine perm_inverse ( n, p ) !*****************************************************************************80 ! !! PERM_INVERSE inverts a permutation "in place". ! ! Modified: ! ! 25 July 2000 ! ! 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 ), parameter :: i4_1 = 1 integer ( kind = 4 ) is integer ( kind = 4 ) p(n) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PERM_INVERSE - 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 ( i4_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 points_delaunay_naive_2d ( node_num, node_xy, maxtri, & triangle_num, triangle_node ) !*****************************************************************************80 ! !! POINTS_DELAUNAY_NAIVE_2D is a naive Delaunay triangulation scheme. ! ! Discussion: ! ! This routine is only suitable as a demonstration code for small ! problems. Its running time is of order NODE_NUM**4. Much faster ! algorithms are available. ! ! Given a set of nodes in the plane, a triangulation is set of ! triples of distinct nodes, forming triangles, so that every ! point within the convex hull of the set of nodes is either ! one of the nodes, or lies on an edge of one or more triangles, ! or lies within exactly one triangle. ! ! A Delaunay triangulation is a triangulation with additional ! properties. ! ! NODE_NUM must be at least 3. ! ! Modified: ! ! 08 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Joseph ORourke, ! Computational Geometry, ! Cambridge University Press, ! Second Edition, 1998, page 187. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. ! ! Input, integer ( kind = 4 ) MAXTRI, the maximum number of triangles. ! ! Output, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles in ! the triangulation. ! ! Output, integer ( kind = 4 ) TRIANGLE_NODE(3,MAXTRI), the indices of ! the triangle nodes. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) maxtri integer ( kind = 4 ) node_num logical flag integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) m real ( kind = 8 ) node_xy(dim_num,node_num) integer ( kind = 4 ) triangle_node(3,maxtri) integer ( kind = 4 ) triangle_num real ( kind = 8 ) xn real ( kind = 8 ) yn real ( kind = 8 ) z(node_num) real ( kind = 8 ) zn triangle_num = 0 if ( node_num < 3 ) then return end if ! ! Compute Z = X*X + Y*Y. ! z(1:node_num) = node_xy(1,1:node_num)**2 + node_xy(2,1:node_num)**2 ! ! For each triple (I,J,K): ! do i = 1, node_num - 2 do j = i+1, node_num do k = i+1, node_num if ( j /= k ) then xn = ( node_xy(2,j) - node_xy(2,i) ) * ( z(k) - z(i) ) & - ( node_xy(2,k) - node_xy(2,i) ) * ( z(j) - z(i) ) yn = ( node_xy(1,k) - node_xy(1,i) ) * ( z(j) - z(i) ) & - ( node_xy(1,j) - node_xy(1,i) ) * ( z(k) - z(i) ) zn = ( node_xy(1,j) - node_xy(1,i) ) & * ( node_xy(2,k) - node_xy(2,i) ) & - ( node_xy(1,k) - node_xy(1,i) ) & * ( node_xy(2,j) - node_xy(2,i) ) flag = ( zn < 0.0D+00 ) if ( flag ) then do m = 1, node_num flag = flag .and. & ( ( node_xy(1,m) - node_xy(1,i) ) * xn & + ( node_xy(2,m) - node_xy(2,i) ) * yn & + ( z(m) - z(i) ) * zn <= 0.0D+00 ) end do end if if ( flag ) then if ( triangle_num < maxtri ) then triangle_num = triangle_num + 1 triangle_node(1:3,triangle_num) = (/ i, j, k /) end if end if end if end do end do end do return end subroutine points_hull_2d ( node_num, node_xy, nval, ival ) !*****************************************************************************80 ! !! POINTS_HULL_2D computes the convex hull of a set of points in 2D. ! ! Modified: ! ! 29 January 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. ! ! Output, integer ( kind = 4 ) NVAL, the number of nodes that lie on the ! convex hull. ! ! Output, integer ( kind = 4 ) IVAL(NODE_NUM). The first NVAL entries of ! IVAL contain the indices of the nodes that form the convex hull, in order. ! These indices are 1-based, not 0-based! ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) ai real ( kind = 8 ) angle_rad_2d real ( kind = 8 ) angmax real ( kind = 8 ) di real ( kind = 8 ) dr integer ( kind = 4 ) i integer ( kind = 4 ) iq integer ( kind = 4 ) ir integer ( kind = 4 ) istart integer ( kind = 4 ) ival(node_num) real ( kind = 8 ) node_xy(dim_num,node_num) integer ( kind = 4 ) nval real ( kind = 8 ) pp(dim_num) if ( node_num < 1 ) then nval = 0 return end if ! ! If NODE_NUM = 1, the hull is the point. ! if ( node_num == 1 ) then nval = 1 ival(1) = 1 return end if ! ! If NODE_NUM = 2, then the convex hull is either the two distinct points, ! or possibly a single (repeated) point. ! if ( node_num == 2 ) then nval = 1 ival(1) = 1 if ( any ( node_xy(1:dim_num,1) /= node_xy(1:dim_num,2) ) ) then nval = nval + 1 ival(2) = 2 end if return end if ! ! Find the leftmost point, and take the bottom-most in a tie. ! Call it "Q". ! iq = 1 do i = 2, node_num if ( node_xy(1,i) < node_xy(1,iq) .or. & ( node_xy(1,i) == node_xy(1,iq) .and. & node_xy(2,i) < node_xy(2,iq) ) ) then iq = i end if end do ! ! Remember the starting point. ! istart = iq nval = 1 ival(1) = iq ! ! For the first point, make a dummy previous point, 1 unit south, ! and call it "P". ! pp(1:2) = (/ node_xy(1,iq), node_xy(2,iq) - 1.0D+00 /) ! ! Now, having old point P, and current point Q, find the new point R ! so the angle PQR is maximal. ! ! Watch out for the possibility that the two nodes are identical. ! do ir = 0 angmax = 0.0D+00 do i = 1, node_num if ( i /= iq .and. & ( node_xy(1,i) /= node_xy(1,iq) .or. & node_xy(2,i) /= node_xy(2,iq) ) ) then ai = angle_rad_2d ( node_xy(1:dim_num,i), node_xy(1:dim_num,iq), pp ) if ( ir == 0 .or. angmax < ai ) then ir = i angmax = ai ! ! In case of ties, choose the nearer point. ! else if ( ir /= 0 .and. ai == angmax ) then di = sqrt ( sum ( & ( node_xy(1:dim_num,i) - node_xy(1:dim_num,iq) )**2 ) ) dr = sqrt ( sum ( & ( node_xy(1:dim_num,ir) - node_xy(1:dim_num,iq) )**2 ) ) if ( di < dr ) then ir = i angmax = ai end if end if end if end do ! ! If we've returned to our starting point, exit. ! if ( ir == istart ) then exit end if if ( node_num < nval + 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_HULL_2D - Fatal error!' write ( *, '(a)' ) ' The algorithm failed.' write ( *, '(a)' ) ' We have not returned to the starting point,' write ( *, '(a)' ) ' and we have already used all the points.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' NODE_NUM = NVAL = ', nval stop end if ! ! Set Q := P, P := R, and repeat. ! nval = nval + 1 ival(nval) = ir pp(1:dim_num) = node_xy(1:dim_num,iq) iq = ir end do return end subroutine points_point_near_naive_nd ( dim_num, nset, pset, p, i_min, d_min ) !*****************************************************************************80 ! !! POINTS_POINT_NEAR_NAIVE_ND finds the nearest point to a given point in ND. ! ! Discussion: ! ! A naive algorithm is used. The distance to every point is calculated, ! in order to determine the smallest. ! ! Modified: ! ! 06 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension. ! ! Input, integer ( kind = 4 ) NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(DIM_NUM,NSET), the points in the set. ! ! Input, real ( kind = 8 ) P(DIM_NUM), the point whose nearest neighbor ! is sought. ! ! Output, integer ( kind = 4 ) I_MIN, the index of the nearest point in ! PSET to P. ! ! Output, real ( kind = 8 ) D_MIN, the distance between P(*) ! and PSET(*,I_MIN). ! implicit none integer ( kind = 4 ) dim_num integer ( kind = 4 ) nset real ( kind = 8 ) d real ( kind = 8 ) d_min integer ( kind = 4 ) i integer ( kind = 4 ) i_min real ( kind = 8 ) p(dim_num) real ( kind = 8 ) pset(dim_num,nset) d_min = huge ( d_min ) i_min = -1 do i = 1, nset d = sum ( ( p(1:dim_num) - pset(1:dim_num,i) )**2 ) if ( d < d_min ) then d_min = d i_min = i end if end do d_min = sqrt ( d_min ) return end function q_measure ( n, z, triangle_order, triangle_num, triangle_node ) !****************************************************************************** ! !! Q_MEASURE determines the triangulated pointset quality measure Q. ! ! Discussion: ! ! The Q measure evaluates the uniformity of the shapes of the triangles ! defined by a triangulated pointset. ! ! For a single triangle T, the value of Q(T) is defined as follows: ! ! TAU_IN = radius of the inscribed circle, ! TAU_OUT = radius of the circumscribed circle, ! ! Q(T) = 2 * TAU_IN / TAU_OUT ! = ( B + C - A ) * ( C + A - B ) * ( A + B - C ) / ( A * B * C ) ! ! where A, B and C are the lengths of the sides of the triangle T. ! ! The Q measure computes the value of Q(T) for every triangle T in the ! triangulation, and then computes the minimum of this ! set of values: ! ! Q_MEASURE = min ( all T in triangulation ) Q(T) ! ! In an ideally regular mesh, all triangles would have the same ! equilateral shape, for which Q = 1. A good mesh would have ! 0.5 < Q. ! ! Given the 2D coordinates of a set of N nodes, stored as Z(1:2,1:N), ! a triangulation is a list of TRIANGLE_NUM triples of node indices that form ! triangles. Generally, a maximal triangulation is expected, namely, ! a triangulation whose image is a planar graph, but for which the ! addition of any new triangle would mean the graph was no longer planar. ! A Delaunay triangulation is a maximal triangulation which maximizes ! the minimum angle that occurs in any triangle. ! ! The code has been modified to 'allow' 6-node triangulations. ! However, no effort is made to actually process the midside nodes. ! Only information from the vertices is used. ! ! Modified: ! ! 07 November 2005 ! ! Author: ! ! Max Gunzburger ! John Burkardt ! ! Reference: ! ! Max Gunzburger, John Burkardt, ! Uniformity Measures for Point Samples in Hypercubes. ! ! Per-Olof Persson, Gilbert Strang, ! A Simple Mesh Generator in MATLAB, ! SIAM Review, ! Volume 46, Number 2, June 2004, pages 329-345. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) Z(2,N), the points. ! ! Input, integer ( kind = 4 ) TRIANGLE_ORDER, the order of the triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM), ! the triangulation. ! ! Output, real ( kind = 8 ) Q_MEASURE, the Q quality measure. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) triangle_num integer ( kind = 4 ) triangle_order integer ( kind = 4 ) a_index real ( kind = 8 ) ab_length integer ( kind = 4 ) b_index real ( kind = 8 ) bc_length integer ( kind = 4 ) c_index real ( kind = 8 ) ca_length real ( kind = 8 ) q real ( kind = 8 ) q_min real ( kind = 8 ) q_measure integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) real ( kind = 8 ) z(2,n) if ( triangle_num < 1 ) then q_measure = -1.0D+00 return end if q_min = huge ( q_min ) do triangle = 1, triangle_num a_index = triangle_node(1,triangle) b_index = triangle_node(2,triangle) c_index = triangle_node(3,triangle) ab_length = sqrt ( & ( z(1,a_index) - z(1,b_index) )**2 & + ( z(2,a_index) - z(2,b_index) )**2 ) bc_length = sqrt ( & ( z(1,b_index) - z(1,c_index) )**2 & + ( z(2,b_index) - z(2,c_index) )**2 ) ca_length = sqrt ( & ( z(1,c_index) - z(1,a_index) )**2 & + ( z(2,c_index) - z(2,a_index) )**2 ) q = ( bc_length + ca_length - ab_length ) & * ( ca_length + ab_length - bc_length ) & * ( ab_length + bc_length - ca_length ) & / ( ab_length * bc_length * ca_length ) q_min = min ( q_min, q ) end do q_measure = q_min return end function r4_uniform_01 ( seed ) !*****************************************************************************80 ! !! R4_UNIFORM_01 returns a unit pseudorandom R4. ! ! Discussion: ! ! An R4 is a real ( kind = 4 ) value. ! ! This routine implements the recursion ! ! seed = 16807 * seed mod ( 2**31 - 1 ) ! r4_uniform_01 = seed / ( 2**31 - 1 ) ! ! The integer arithmetic never requires more than 32 bits, ! including a sign bit. ! ! If the initial seed is 12345, then the first three computations are ! ! Input Output R4_UNIFORM_01 ! SEED SEED ! ! 12345 207482415 0.096616 ! 207482415 1790989824 0.833995 ! 1790989824 2035175616 0.947702 ! ! Modified: ! ! 11 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley Interscience, page 95, 1998. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 4 ) R4_UNIFORM_01, a new pseudorandom variate, ! strictly between 0 and 1. ! implicit none integer ( kind = 4 ) k integer ( kind = 4 ) seed real ( kind = 4 ) r4_uniform_01 if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R4_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if r4_uniform_01 = real ( seed, kind = 4 ) * 4.656612875E-10 return end function r8_uniform_01 ( seed ) !*****************************************************************************80 ! !! R8_UNIFORM_01 returns a unit pseudorandom R8. ! ! Discussion: ! ! An R8 is a real ( kind = 8 ) value. ! ! For now, the input quantity SEED is an integer ( kind = 4 ) variable. ! ! This routine implements the recursion ! ! seed = 16807 * seed mod ( 2**31 - 1 ) ! r8_uniform_01 = seed / ( 2**31 - 1 ) ! ! The integer arithmetic never requires more than 32 bits, ! including a sign bit. ! ! If the initial seed is 12345, then the first three computations are ! ! Input Output R8_UNIFORM_01 ! SEED SEED ! ! 12345 207482415 0.096616 ! 207482415 1790989824 0.833995 ! 1790989824 2035175616 0.947702 ! ! Modified: ! ! 05 July 2006 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley Interscience, page 95, 1998. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R8_UNIFORM_01, a new pseudorandom variate, ! strictly between 0 and 1. ! implicit none integer ( kind = 4 ) k real ( kind = 8 ) r8_uniform_01 integer ( kind = 4 ) seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if ! ! Although SEED can be represented exactly as a 32 bit integer, ! it generally cannot be represented exactly as a 32 bit real number! ! r8_uniform_01 = real ( seed, kind = 8 ) * 4.656612875D-10 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 ). ! ! Modified: ! ! 13 March 2005 ! ! 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 integer ( kind = 4 ), parameter :: ndim = 2 real ( kind = 8 ) a(ndim,n) real ( kind = 8 ) a_temp(ndim) integer ( kind = 4 ) ierror integer ( kind = 4 ) iget integer ( kind = 4 ) iput integer ( kind = 4 ) istart integer ( kind = 4 ) p(n) call perm_check ( n, p, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R82VEC_PERMUTE - Fatal error!' write ( *, '(a)' ) ' The input array does not represent' write ( *, '(a)' ) ' a proper permutation. In particular, the' write ( *, '(a,i8)' ) ' array is missing the value ', ierror stop end if ! ! 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:ndim) = a(1:ndim,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)' ) 'R822VEC_PERMUTE - Fatal error!' write ( *, '(a)' ) ' A permutation index is out of range.' write ( *, '(a,i8,a,i8)' ) ' P(', iput, ') = ', iget stop end if if ( iget == istart ) then a(1:ndim,iput) = a_temp(1:ndim) exit end if a(1:ndim,iput) = a(1:ndim,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. ! ! 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 r8mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_PRINT prints an R8MAT. ! ! Modified: ! ! 20 May 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows in A. ! ! Input, integer ( kind = 4 ) N, the number of columns in A. ! ! Input, real ( kind = 8 ) A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title to be printed. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) integer ( kind = 4 ), parameter :: i4_1 = 1 character ( len = * ) title call r8mat_print_some ( m, n, a, i4_1, i4_1, m, n, title ) return end subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_PRINT_SOME prints some of an R8MAT. ! ! Modified: ! ! 26 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i7,7x)') j end do write ( *, '('' Col '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 if ( a(i,j) == real ( int ( a(i,j) ), kind = 8 ) ) then write ( ctemp(j2), '(f8.0,6x)' ) a(i,j) else write ( ctemp(j2), '(g14.6)' ) a(i,j) end if end do write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j), j = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine r8mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) integer ( kind = 4 ), parameter :: i4_1 = 1 character ( len = * ) title call r8mat_transpose_print_some ( m, n, a, i4_1, i4_1, m, n, title ) return end subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2 integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7,7x)') i end do write ( *, '('' Row '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,1x,5a14)' ) j, ( ctemp(i), i = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine r8tris2 ( node_num, node_xy, triangle_num, triangle_node, & triangle_neighbor ) !*****************************************************************************80 ! !! R8TRIS2 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. ! ! Modified: ! ! 25 August 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! 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 ) NODE_NUM, the number of nodes. ! ! Input/output, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates ! of the nodes. On output, the vertices have been sorted into ! dictionary order. ! ! Output, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles in the ! triangulation; TRIANGLE_NUM is equal to 2*NODE_NUM - NB - 2, where NB is ! the number of boundary vertices. ! ! Output, integer ( kind = 4 ) TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes that ! make up each triangle. The elements are indices of P. The vertices of ! the triangles are in counter clockwise order. ! ! Output, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_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; ! TRIANGLE_NEIGHBOR(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 ) node_num real ( kind = 8 ) cmax integer ( kind = 4 ) e integer ( kind = 4 ) i integer ( kind = 4 ), parameter :: i4_3 = 3 integer ( kind = 4 ) ierr integer ( kind = 4 ) indx(node_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 ) node_xy(dim_num,node_num) integer ( kind = 4 ) redg integer ( kind = 4 ) rtri integer ( kind = 4 ) stack(node_num) integer ( kind = 4 ) t real ( kind = 8 ) tol integer ( kind = 4 ) top integer ( kind = 4 ) triangle_neighbor(3,node_num*2) integer ( kind = 4 ) triangle_num integer ( kind = 4 ) triangle_node(3,node_num*2) tol = 100.0D+00 * epsilon ( tol ) ierr = 0 ! ! Sort the vertices by increasing (x,y). ! call r82vec_sort_heap_index_a ( node_num, node_xy, indx ) call r82vec_permute ( node_num, node_xy, indx ) ! ! Make sure that the data nodes are "reasonably" distinct. ! m1 = 1 do i = 2, node_num m = m1 m1 = i k = 0 do j = 1, dim_num cmax = max ( abs ( node_xy(j,m) ), abs ( node_xy(j,m1) ) ) if ( tol * ( cmax + 1.0D+00 ) & < abs ( node_xy(j,m) - node_xy(j,m1) ) ) then k = j exit end if end do if ( k == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8TRIS2 - 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)' ) ' NODE_XY(M) = ', node_xy(1:dim_num,m) write ( *, '(a,2g14.6)' ) ' NODE_XY(M1) = ', node_xy(1:dim_num,m1) ierr = 224 stop end if end do ! ! Starting from nodes 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 ( node_num < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8TRIS2 - Fatal error!' ierr = 225 stop end if m = j lr = lrline ( node_xy(1,m), node_xy(2,m), node_xy(1,m1), & node_xy(2,m1), node_xy(1,m2), node_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. ! triangle_num = j - 2 if ( lr == -1 ) then triangle_node(1,1) = m1 triangle_node(2,1) = m2 triangle_node(3,1) = m triangle_neighbor(3,1) = -3 do i = 2, triangle_num m1 = m2 m2 = i+1 triangle_node(1,i) = m1 triangle_node(2,i) = m2 triangle_node(3,i) = m triangle_neighbor(1,i-1) = -3 * i triangle_neighbor(2,i-1) = i triangle_neighbor(3,i) = i - 1 end do triangle_neighbor(1,triangle_num) = -3 * triangle_num - 1 triangle_neighbor(2,triangle_num) = -5 ledg = 2 ltri = triangle_num else triangle_node(1,1) = m2 triangle_node(2,1) = m1 triangle_node(3,1) = m triangle_neighbor(1,1) = -4 do i = 2, triangle_num m1 = m2 m2 = i+1 triangle_node(1,i) = m2 triangle_node(2,i) = m1 triangle_node(3,i) = m triangle_neighbor(3,i-1) = i triangle_neighbor(1,i) = -3 * i - 3 triangle_neighbor(2,i) = i - 1 end do triangle_neighbor(3,triangle_num) = -3 * triangle_num triangle_neighbor(2,1) = -3 * triangle_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, node_num m = i m1 = triangle_node(ledg,ltri) if ( ledg <= 2 ) then m2 = triangle_node(ledg+1,ltri) else m2 = triangle_node(1,ltri) end if lr = lrline ( node_xy(1,m), node_xy(2,m), node_xy(1,m1), & node_xy(2,m1), node_xy(1,m2), node_xy(2,m2), 0.0D+00 ) if ( 0 < lr ) then rtri = ltri redg = ledg ltri = 0 else l = -triangle_neighbor(ledg,ltri) rtri = l / 3 redg = mod ( l, i4_3 ) + 1 end if call vbedg ( node_xy(1,m), node_xy(2,m), node_num, node_xy, triangle_num, & triangle_node, triangle_neighbor, ltri, ledg, rtri, redg ) n = triangle_num + 1 l = -triangle_neighbor(ledg,ltri) do t = l / 3 e = mod ( l, i4_3 ) + 1 l = -triangle_neighbor(e,t) m2 = triangle_node(e,t) if ( e <= 2 ) then m1 = triangle_node(e+1,t) else m1 = triangle_node(1,t) end if triangle_num = triangle_num + 1 triangle_neighbor(e,t) = triangle_num triangle_node(1,triangle_num) = m1 triangle_node(2,triangle_num) = m2 triangle_node(3,triangle_num) = m triangle_neighbor(1,triangle_num) = t triangle_neighbor(2,triangle_num) = triangle_num - 1 triangle_neighbor(3,triangle_num) = triangle_num + 1 top = top + 1 if ( node_num < top ) then ierr = 8 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8TRIS2 - Fatal error!' write ( *, '(a)' ) ' Stack overflow.' stop end if stack(top) = triangle_num if ( t == rtri .and. e == redg ) then exit end if end do triangle_neighbor(ledg,ltri) = -3 * n - 1 triangle_neighbor(2,n) = -3 * triangle_num - 2 triangle_neighbor(3,triangle_num) = -l ltri = n ledg = 2 call swapec ( m, top, ltri, ledg, node_num, node_xy, triangle_num, & triangle_node, triangle_neighbor, stack, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8TRIS2 - Fatal error!' write ( *, '(a)' ) ' Error return from SWAPEC.' stop end if end do ! ! Now account for the sorting that we did. ! do i = 1, 3 do j = 1, triangle_num triangle_node(i,j) = indx ( triangle_node(i,j) ) end do end do call perm_inverse ( node_num, indx ) call r82vec_permute ( node_num, node_xy, indx ) return end subroutine r8vec_bracket ( n, x, xval, left, right ) !*****************************************************************************80 ! !! R8VEC_BRACKET searches a sorted R8VEC for successive brackets of a value. ! ! Discussion: ! ! If the values in the vector are thought of as defining intervals ! on the real line, then this routine searches for the interval ! nearest to or containing the given value. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, length of input array. ! ! Input, real ( kind = 8 ) X(N), an array that has been sorted into ! ascending order. ! ! Input, real ( kind = 8 ) XVAL, a value to be bracketed. ! ! Output, integer ( kind = 4 ) LEFT, RIGHT, the results of the search. ! Either: ! XVAL < X(1), when LEFT = 1, RIGHT = 2; ! X(N) < XVAL, when LEFT = N-1, RIGHT = N; ! or ! X(LEFT) <= XVAL <= X(RIGHT). ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) left integer ( kind = 4 ) right real ( kind = 8 ) x(n) real ( kind = 8 ) xval do i = 2, n - 1 if ( xval < x(i) ) then left = i - 1 right = i return end if end do left = n - 1 right = n return end subroutine sort_heap_external ( n, indx, i, j, isgn ) !*****************************************************************************80 ! !! SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. ! ! Discussion: ! ! The actual list of data is not passed to the routine. Hence this ! routine may be used to sort integers, real ( kind = 8 )s, numbers, names, ! dates, shoe sizes, and so on. After each call, the routine asks ! the user to compare or interchange two items, until a special ! return value signals that the sorting is completed. ! ! Modified: ! ! 05 February 2004 ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of items to be sorted. ! ! Input/output, integer ( kind = 4 ) INDX, the main communication signal. ! ! The user must set INDX to 0 before the first call. ! Thereafter, the user should not change the value of INDX until ! the sorting is done. ! ! On return, if INDX is ! ! greater than 0, ! * interchange items I and J; ! * call again. ! ! less than 0, ! * compare items I and J; ! * set ISGN = -1 if I < J, ISGN = +1 if J < I; ! * call again. ! ! equal to 0, the sorting is done. ! ! Output, integer ( kind = 4 ) I, J, the indices of two items. ! On return with INDX positive, elements I and J should be interchanged. ! On return with INDX negative, elements I and J should be compared, and ! the result reported in ISGN on the next call. ! ! Input, integer ( kind = 4 ) ISGN, results of comparison of elements I ! and J. (Used only when the previous call returned INDX less than 0). ! ISGN <= 0 means I is less than or equal to J; ! 0 <= ISGN means I is greater than or equal to J. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ), save :: i_save = 0 integer ( kind = 4 ) indx integer ( kind = 4 ) isgn integer ( kind = 4 ) j integer ( kind = 4 ), save :: j_save = 0 integer ( kind = 4 ), save :: k = 0 integer ( kind = 4 ), save :: k1 = 0 integer ( kind = 4 ) n integer ( kind = 4 ), save :: n1 = 0 ! ! INDX = 0: This is the first call. ! if ( indx == 0 ) then i_save = 0 j_save = 0 k = n / 2 k1 = k n1 = n ! ! INDX < 0: The user is returning the results of a comparison. ! else if ( indx < 0 ) then if ( indx == -2 ) then if ( isgn < 0 ) then i_save = i_save + 1 end if j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return end if if ( 0 < isgn ) then indx = 2 i = i_save j = j_save return end if if ( k <= 1 ) then if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 end if i = i_save j = j_save return end if k = k - 1 k1 = k ! ! 0 < INDX, the user was asked to make an interchange. ! else if ( indx == 1 ) then k1 = k end if do i_save = 2 * k1 if ( i_save == n1 ) then j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return else if ( i_save <= n1 ) then j_save = i_save + 1 indx = -2 i = i_save j = j_save return end if if ( k <= 1 ) then exit end if k = k - 1 k1 = k end do if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 i = i_save j = j_save else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 i = i_save j = j_save end if return end subroutine swapec ( i, top, btri, bedg, node_num, node_xy, triangle_num, & triangle_node, triangle_neighbor, 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. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! 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 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input/output, integer ( kind = 4 ) TRIANGLE_NODE(3,TRIANGLE_NUM), the ! triangle incidence list. May be updated on output because of swaps. ! ! Input/output, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_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 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 ( kind = 4 ) IERR is set to 8 for abnormal return. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ), parameter :: i4_1 = 1 integer ( kind = 4 ), parameter :: i4_3 = 3 integer ( kind = 4 ) node_num integer ( kind = 4 ) triangle_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 ) node_xy(dim_num,node_num) integer ( kind = 4 ) r integer ( kind = 4 ) s integer ( kind = 4 ) stack(node_num) integer ( kind = 4 ) swap integer ( kind = 4 ) t integer ( kind = 4 ) top integer ( kind = 4 ) triangle_node(3,triangle_num) integer ( kind = 4 ) triangle_neighbor(3,triangle_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 = node_xy(1,i) y = node_xy(2,i) do if ( top <= 0 ) then exit end if t = stack(top) top = top - 1 if ( triangle_node(1,t) == i ) then e = 2 b = triangle_node(3,t) else if ( triangle_node(2,t) == i ) then e = 3 b = triangle_node(1,t) else e = 1 b = triangle_node(2,t) end if a = triangle_node(e,t) u = triangle_neighbor(e,t) if ( triangle_neighbor(1,u) == t ) then f = 1 c = triangle_node(3,u) else if ( triangle_neighbor(2,u) == t ) then f = 2 c = triangle_node(1,u) else f = 3 c = triangle_node(2,u) end if swap = diaedg ( x, y, node_xy(1,a), node_xy(2,a), node_xy(1,c), & node_xy(2,c), node_xy(1,b), node_xy(2,b) ) if ( swap == 1 ) then em1 = i4_wrap ( e - i4_1, i4_1, i4_3 ) ep1 = i4_wrap ( e + i4_1, i4_1, i4_3 ) fm1 = i4_wrap ( f - i4_1, i4_1, i4_3 ) fp1 = i4_wrap ( f + i4_1, i4_1, i4_3 ) triangle_node(ep1,t) = c triangle_node(fp1,u) = i r = triangle_neighbor(ep1,t) s = triangle_neighbor(fp1,u) triangle_neighbor(ep1,t) = u triangle_neighbor(fp1,u) = t triangle_neighbor(e,t) = s triangle_neighbor(f,u) = r if ( 0 < triangle_neighbor(fm1,u) ) then top = top + 1 stack(top) = u end if if ( 0 < s ) then if ( triangle_neighbor(1,s) == u ) then triangle_neighbor(1,s) = t else if ( triangle_neighbor(2,s) == u ) then triangle_neighbor(2,s) = t else triangle_neighbor(3,s) = t end if top = top + 1 if ( node_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 < triangle_neighbor(ee,tt) ) tt = triangle_neighbor(ee,tt) if ( triangle_node(1,tt) == a ) then ee = 3 else if ( triangle_node(2,tt) == a ) then ee = 1 else ee = 2 end if end do triangle_neighbor(ee,tt) = l end if if ( 0 < r ) then if ( triangle_neighbor(1,r) == t ) then triangle_neighbor(1,r) = u else if ( triangle_neighbor(2,r) == t ) then triangle_neighbor(2,r) = u else triangle_neighbor(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 < triangle_neighbor(ee,tt) ) tt = triangle_neighbor(ee,tt) if ( triangle_node(1,tt) == b ) then ee = 3 else if ( triangle_node(2,tt) == b ) then ee = 1 else ee = 2 end if end do triangle_neighbor(ee,tt) = l end if end if end do return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 26 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer ( kind = 4 ) d integer ( kind = 4 ) h integer ( kind = 4 ) m integer ( kind = 4 ) mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer ( kind = 4 ) n integer ( kind = 4 ) s integer ( kind = 4 ) values(8) integer ( kind = 4 ) y call date_and_