cvt_mesh.f90

一个三维有限元网格自动剖分源代码程序

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 ( )

  stop
end
subroutine 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 )

  return
end
subroutine 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 )

  return
end
subroutine 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) )