program main!*****************************************************************************80!!! MAIN is the main program for CVT_TET_MESH.!!  Discussion:!!    CVT_TET_MESH uses CVT sampling on problems from TEST_TET_MESH.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    28 October 2005!!  Author:!!    John Burkardt!  implicit none  integer test  integer test_num  call timestamp ( )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) 'CVT_TET_MESH:'  write ( *, '(a)' ) '  FORTRAN90 version'  write ( *, '(a)' ) '  Apply simple CVT sampling routines to produce'  write ( *, '(a)' ) '  a set of sample points in regions from'  write ( *, '(a)' ) '  the TEST_TET_MESH package.'!!  How many test cases are available?!  call p00_test_num ( test_num )  write ( *, '(a)' ) ' '  write ( *, '(a,i6)' ) '  Number of test cases = ', test_num!!  Run test 1 on the test cases.  This simply computes a CVT in the!  interior of the test region.!  if ( .false. ) then    do test = 1, test_num      call test01 ( test )    end do  end if!!  Run test 2 on the test cases.  This computes a nonstandard CVT in the!  interior and on the boundary of the test region.!  if ( .false. ) then    do test = 1, test_num      call test02 ( test )    end do  end if!!  Run test 3 on the test cases.  This computes a nonstandard CVT in the!  interior and on the boundary of the test region, plus fixed points.!  if ( .true. ) then    do test = 2, test_num      call test03 ( test )    end do  end if  write ( *, '(a)' ) ' '  write ( *, '(a)' ) 'CVT_TET_MESH:'  write ( *, '(a)' ) '  Normal end of execution.'  write ( *, '(a)' ) ' '  call timestamp ( )  stopendsubroutine test01 ( test )!*****************************************************************************80!!! TEST01 creates a standard CVT in a given test region.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    06 August 2005!!  Author:!!    John Burkardt!!  Parameters:!!    Input, integer TEST, the index of the problem to be treated.!  implicit none  integer, parameter :: dim_num = 3  logical :: comment = .false.  integer, allocatable, dimension ( : ) :: count  real ( kind = 8 ) energy  character ( len = 80 ) file_eps_name  character ( len = 80 ) file_txt_name  integer iteration  integer, parameter :: iteration_max = 40  integer j  integer, parameter :: n = 400  integer, allocatable, dimension ( : ) :: nearest  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point_new  real ( kind = 8 ), allocatable, dimension ( :, : ) :: sample  integer, parameter :: sample_num = 100000  integer seed  integer :: seed_start = 123456789  integer test  character ( len = 80 ) title  seed = seed_start  energy = -1.0D+00  write ( *, '(a)' ) ' '  write ( *, '(a)' ) 'TEST01'  write ( *, '(a)' ) '  Compute a standard CVT in the interior of'  write ( *, '(a)' ) '  a test region from TEST_TET_MESH.'  call p00_title ( test, title )  write ( *, '(a)' ) ' '  write ( *, '(a, a)' ) '  Title:             "', trim ( title ) // '"'  call p00_header ( test )!!  Initialize the sampling points.!  allocate ( count(1:n) )  allocate ( nearest(1:sample_num) )  allocate ( point(1:dim_num,1:n) )  allocate ( point_new(1:dim_num,1:n) )  allocate ( sample(1:dim_num,1:sample_num) )  call p00_sample ( test, n, seed, point )  call r8mat_transpose_print_some ( dim_num, n, point, &    1, 1, 2, 10, '  Initial points (first 10 only)' )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) '  Estimated Voronoi energy:'  write ( *, '(a)' ) ' '!!  Carry out an iteration!  do iteration = 1, iteration_max    call p00_sample ( test, sample_num, seed, sample )!!  Find the nearest cell generator.!    call find_closest ( dim_num, n, sample_num, point, sample, nearest )!!  Add X to the averaging data for CELL_GENERATOR(*,NEAREST).!    point_new(1:dim_num,1:n) = 0.0D+00    count(1:n) = 0    energy = 0.0D+00    do j = 1, sample_num      energy = energy &        + sum ( ( point(1:dim_num,nearest(j)) - sample(1:dim_num,j) )**2 )      point_new(1:dim_num,nearest(j)) = point_new(1:dim_num,nearest(j)) &        + sample(1:dim_num,j)      count(nearest(j)) = count(nearest(j)) + 1    end do!!  Compute the new generators.!    do j = 1, n      if ( count(j) /= 0 ) then        point_new(1:dim_num,j) = point_new(1:dim_num,j) &          / real ( count(j), kind = 8 )      end if    end do    energy = energy / real ( sample_num, kind = 8 )    write ( *, '(2x,i6,2x,g14.6)' ) iteration, energy!!  Update.!    point(1:dim_num,1:n) = point_new(1:dim_num,1:n)  end do  if ( test < 10 ) then    write ( file_txt_name, '(a,i1,a)' ) 'cvt_p0', test, '_test01.txt'  else    write ( file_txt_name, '(a,i2,a)' ) 'cvt_p', test, '_test01.txt'  end if  call cvt_write ( dim_num, n, seed_start, seed, sample_num, iteration_max, &    energy, point, file_txt_name, comment )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) '  TEST01: wrote data to "' &    // trim ( file_txt_name ) // '".'  deallocate ( count )  deallocate ( nearest )  deallocate ( point )  deallocate ( point_new )   deallocate ( sample )  returnendsubroutine test02 ( test )!*****************************************************************************80!!! TEST02 creates a CVT dataset, modified so some points are on the boundary.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    06 August 2005!!  Author:!!    John Burkardt!!  Parameters:!!    Input, integer TEST, the index of the problem to be treated.!  implicit none  integer, parameter :: dim_num = 3  logical :: comment = .false.  integer, allocatable, dimension ( : ) :: count  real ( kind = 8 ) energy  character ( len = 80 ) file_txt_name  real ( kind = 8 ) :: h = 0.05D+00  real ( kind = 8 ) hi(dim_num)  integer iteration  integer, parameter :: iteration_max = 40  integer j  real ( kind = 8 ) lo(dim_num)  integer, parameter :: n = 400  integer, allocatable, dimension ( : ) :: nearest  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point_new  real ( kind = 8 ), allocatable, dimension ( :, : ) :: sample  integer, parameter :: sample_num = 100000  real ( kind = 8 ) scale  integer seed  integer :: seed_start = 123456789  integer test  character ( len = 80 ) title  seed = seed_start  energy = -1.0D+00  write ( *, '(a)' ) ' '  write ( *, '(a)' ) 'TEST02'  write ( *, '(a)' ) '  Try to compute a nonstandard CVT in the interior and'  write ( *, '(a)' ) '  on the boundary of a test region from TEST_TET_MESH.'  call p00_title ( test, title )  write ( *, '(a)' ) ' '  write ( *, '(a, a)' ) '  Title:             "', trim ( title ) // '"'  call p00_header ( test )!!  Allocate space for some arrays.!  allocate ( count(1:n) )  allocate ( nearest(1:sample_num) )  allocate ( point(1:dim_num,1:n) )  allocate ( point_new(1:dim_num,1:n) )  allocate ( sample(1:dim_num,1:sample_num) )!!  Determine the amount by which the region should be expanded.!  call p00_box ( test, lo, hi )  scale = maxval ( hi(1:3) - lo(1:3) )  h = 0.05D+00 * scale!!  Initialize the sampling points.!  call p00_sample ( test, n, seed, point )  call r8mat_transpose_print_some ( dim_num, n, point, &    1, 1, 2, 10, '  Initial points (first 10 only)' )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) '  Estimated Voronoi energy:'  write ( *, '(a)' ) ' '!!  Carry out an iteration!  do iteration = 1, iteration_max!!  Get sample points from a "slightly enlarged" region.!    call p00_sample_h1 ( test, sample_num, h, seed, sample )!!  Find the nearest cell generator.!    call find_closest ( dim_num, n, sample_num, point, sample, nearest )!!  Add X to the averaging data for CELL_GENERATOR(*,NEAREST).!    point_new(1:dim_num,1:n) = 0.0D+00    count(1:n) = 0    energy = 0.0D+00    do j = 1, sample_num      energy = energy &        + sum ( ( point(1:dim_num,nearest(j)) - sample(1:dim_num,j) )**2 )      point_new(1:dim_num,nearest(j)) = point_new(1:dim_num,nearest(j)) &        + sample(1:dim_num,j)      count(nearest(j)) = count(nearest(j)) + 1    end do!!  Compute the new generators.!    do j = 1, n      if ( count(j) /= 0 ) then        point_new(1:dim_num,j) = point_new(1:dim_num,j) &          / real ( count(j), kind = 8 )      end if    end do!!  Project generators back into region.!    call p00_boundary_project ( test, n, point_new )    energy = energy / real ( sample_num, kind = 8 )    write ( *, '(2x,i6,2x,g14.6)' ) iteration, energy!!  Update.!    point(1:dim_num,1:n) = point_new(1:dim_num,1:n)  end do  if ( test < 10 ) then    write ( file_txt_name, '(a,i1,a)' ) 'cvt_p0', test, '_test02.txt'  else    write ( file_txt_name, '(a,i2,a)' ) 'cvt_p', test, '_test02.txt'  end if  call cvt_write ( dim_num, n, seed_start, seed, sample_num, iteration_max, &    energy, point, file_txt_name, comment )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) '  TEST02: wrote data to "' &    // trim ( file_txt_name ) // '".'  deallocate ( count )  deallocate ( nearest )  deallocate ( point )  deallocate ( point_new )   deallocate ( sample )  returnendsubroutine test03 ( test )!*****************************************************************************80!!! TEST03 creates a CVT dataset with points on the boundary or set by user.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    26 August 2005!!  Author:!!    John Burkardt!!  Parameters:!!    Input, integer TEST, the index of the problem to be treated.!  implicit none  integer, parameter :: dim_num = 3  real ( kind = 8 ) box_volume  logical :: comment = .false.  integer, allocatable, dimension ( : ) :: count  real ( kind = 8 ) energy  character ( len = 80 ) file_txt_name  real ( kind = 8 ), allocatable, dimension ( :, : ) :: fixed  integer fixed_num  real ( kind = 8 ) h  real ( kind = 8 ) hi(dim_num)  integer iteration  integer, parameter :: iteration_max = 20  integer j  real ( kind = 8 ) lo(dim_num)  integer, parameter :: n = 400  integer, allocatable, dimension ( : ) :: nearest  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point_new  real ( kind = 8 ), allocatable, dimension ( :, : ) :: sample  integer, parameter :: sample_num = 1000000  real ( kind = 8 ) scale  integer seed  integer :: seed_start = 123456789  integer test  character ( len = 80 ) title  call p00_fixed_num ( test, fixed_num )  seed = seed_start  energy = -1.0D+00  write ( *, '(a)' ) ' '  write ( *, '(a)' ) 'TEST03'  write ( *, '(a)' ) '  Try to compute a nonstandard CVT in the interior and'  write ( *, '(a)' ) '  on the boundary of a test region from TEST_TET_MESH.'  call p00_title ( test, title )  write ( *, '(a)' ) ' '  write ( *, '(a, a)' ) '  Title:             "', trim ( title ) // '"'  call p00_header ( test )!!  Allocate space for some arrays.!  allocate ( count(1:n) )