meshless.f90

meshles是不采用网格

program main

!*****************************************************************************80
!
!! MAIN is the main program for the meshless basis function routines.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    14 April 2006
!
!  Author:
!
!    Lili Ju
!
  implicit none

  real ( kind = 8 ), allocatable :: a(:)
  real ( kind = 8 ), dimension ( 2 ) :: alpha = (/ 0.0D+00, 1.0D+00 /)
  real ( kind = 8 ), allocatable :: b(:)
  integer, allocatable :: base(:)
  real ( kind = 8 ), allocatable :: basis_center(:,:)
  character ( len = 80 ) :: basis_file = 'basis.dat'
  integer :: basis_m = 0
  real ( kind = 8 ), allocatable :: basis_radius(:)
  real ( kind = 8 ) :: basis_radius_factor = 1.0D+00
  real ( kind = 8 ), dimension ( 2 ) :: beta = (/ 0.0D+00, 1.0D+00 /)
  integer :: center_choice = 0
  character ( len = 80 ) command
  real ( kind = 8 ), allocatable :: cvt_center(:,:)
  character ( len = 80 ) :: cvt_file = 'cvt.dat'
  integer cvt_m
  real ( kind = 8 ), allocatable :: cvt_radius(:)
  real ( kind = 8 ) :: cvt_radius_factor = 1.0D+00
  integer :: density_function = 0
  character ( len = 80 ) :: halton_file = 'halton.dat'
  real ( kind = 8 ), allocatable :: halton_points(:,:)
  integer i_temp
  integer ierror
  integer ios
  integer last
  integer :: maxit = 0
  integer :: myrank = 0
  integer :: n = 0
  integer ncolumn
  integer :: m = 2
  integer :: ns = 0
  integer :: radius_choice = 0
  logical s_eqi
  integer :: skip = 0
  character ( len = 80 ) :: uniform_file = 'uniform.dat'
  real ( kind = 8 ), allocatable :: uniform_points(:,:)
  character ( len = 80 ) what

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'MESHLESS'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Meshless basis function generation.'
!
!  Initialize the random number generator.
!
  call set_random_seed ( myrank )
!
!  Read a command, do a command.
!
  do

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Enter command (or HELP for a list)'

    do

      read ( *, '(a)', iostat = ios ) command
      write ( *, '(a)' ) '  Command: " ' // trim ( command ) // '".'

      if ( command(1:1) /= '#' ) then
        exit
      end if

    end do

    call s_blank_delete ( command )

    if ( ios /= 0 ) then
      exit
    end if

    if ( s_eqi ( command, 'BASIS_MAKE' ) ) then

      if ( m <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'BASIS_MAKE - Error!'
        write ( *, '(a)' ) '  The spatial dimension must be positive!'
        write ( *, '(a)' ) '  Enter the spatial dimension M with the command'
        write ( *, '(a)' ) '    "M = ???"'
        cycle
      end if

      if ( n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'BASIS_MAKE - Error!'
        write ( *, '(a)' ) '  The number of points N must be positive!'
        write ( *, '(a)' ) '  Enter the number of points with the command'
        write ( *, '(a)' ) '    "N = ???"'
        cycle
      end if

      basis_m = 7 * int ( sqrt ( real ( n, kind = 8 ) ) )

      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 'Mysterious parameter BASIS_M = ', basis_m

      do

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Initial center point determination:'
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  1: use a uniform distribution;'
        write ( *, '(a)' ) '  2: use a Halton distribution.'

        read ( *, * ) center_choice

        if ( center_choice == 1 ) then

          call uniform_make ( m, n, a, b, density_function, &
            basis_center )
          exit

        else if ( center_choice == 2 ) then

          call halton_generate ( m, n, skip, base, basis_center )
          exit
        end if

      end do

      do

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Initial radius determination:'
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  1: radius method 1;'
        write ( *, '(a)' ) '  2: radius method 2.'

        read ( *, * ) radius_choice

        if ( radius_choice == 1 ) then

          if ( m <= 0 ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'RADIUS_MAKE_1 - Error!'
            write ( *, '(a)' ) '  The spatial dimension must be positive!'
            write ( *, '(a)' ) &
              '  Enter the spatial dimension M with the command'
            write ( *, '(a)' ) '    "M = ???"'
            exit
          end if

          call radius_make_1 ( m, a, b, n, basis_m, density_function, &
            basis_center, basis_radius )

          exit

        else if ( radius_choice == 2 ) then

          if ( m <= 0 ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'RADIUS_MAKE_2 - Error!'
            write ( *, '(a)' ) '  The spatial dimension must be positive!'
            write ( *, '(a)' ) &
              '  Enter the spatial dimension M with the command'
            write ( *, '(a)' ) '    "M = ???"'
            exit
          end if

          call radius_make_2 ( m, a, b, n, basis_m, basis_center, &
            basis_radius )

          exit
        end if

      end do

    else if ( s_eqi ( command, 'BASIS_OVERLAP' ) ) then

      call basis_overlap ( m, n, basis_center, basis_radius )

    else if ( s_eqi ( command, 'BASIS_PLOT' ) ) then

      call basis_plot ( m, n, basis_center, basis_radius )

    else if ( s_eqi ( command, 'BASIS_READ' ) ) then

      call file_column_count ( basis_file, ncolumn )

      m = ncolumn - 1

      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 'Setting M to ', m

      call file_line_count ( basis_file, n )

      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 'Setting number of center points to ', n

      call basis_read ( m, n, basis_center, basis_radius, basis_file )

      call basis_box ( m, n, basis_center, basis_radius, a, b )

    else if ( s_eqi ( command, 'BASIS_SCALE' ) ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Enter a scale factor to apply to the radius'
      write ( *, '(a)' ) 'of the current basis function set:'

      read ( *, *, iostat = ios ) basis_radius_factor
      if ( ios /= 0 ) then
        exit
      end if

      basis_radius(1:n) = basis_radius_factor * basis_radius(1:n)

    else if ( s_eqi ( command, 'BASIS_WRITE' ) ) then

      if ( n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Use the command radius_make first!'
        cycle
      end if

      call basis_write ( m, n, basis_center, basis_radius, basis_file )

    else if ( s_eqi ( command, 'CENTER_PLOT' ) ) then

      call center_plot ( m, n, basis_center, a, b )

    else if ( s_eqi ( command, 'CVT_MAKE' ) ) then

      if ( m <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CVT_MAKE - Error!'
        write ( *, '(a)' ) '  The spatial dimension must be positive!'
        write ( *, '(a)' ) '  Enter the spatial dimension M with the command'
        write ( *, '(a)' ) '    "M = ???"'
        cycle
      end if

      if ( n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CVT_MAKE - Error!'
        write ( *, '(a)' ) '  The number of points N must be positive!'
        write ( *, '(a)' ) '  Enter the number of points with the command'
        write ( *, '(a)' ) '    "N = ???"'
        cycle
      end if

      cvt_m = 7 * int ( sqrt ( real ( n, kind = 8 ) ) )

      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 'Mysterious parameter CVT_M = ', cvt_m

      do

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Choose a basis center initialization:'
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  1: use a uniform distribution;'
        write ( *, '(a)' ) '  2: use a Halton distribution.'

        read ( *, * ) center_choice

        if ( center_choice == 1 ) then
          call uniform_make ( m, n, a, b, density_function, cvt_center )
          exit
        else if ( center_choice == 2 ) then

          call halton_generate ( m, n, skip, base, cvt_center )
          exit
        end if

      end do

      do

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Initial radius determination:'
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  1: radius method 1;'
        write ( *, '(a)' ) '  2: radius method 2.'

        read ( *, * ) radius_choice

        if ( radius_choice == 1 ) then

          if ( m <= 0 ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'RADIUS_MAKE_2 - Error!'
            write ( *, '(a)' ) '  The spatial dimension must be positive!'
            write ( *, '(a)' ) &
              '  Enter the spatial dimension M with the command'
            write ( *, '(a)' ) '    "M = ???"'
            exit
          end if

          call radius_make_1 ( m, a, b, n, cvt_m, density_function, &
            cvt_center, cvt_radius )
          exit

        else if ( radius_choice == 2 ) then

          if ( m <= 0 ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'RADIUS_MAKE_2 - Error!'
            write ( *, '(a)' ) '  The spatial dimension must be positive!'
            write ( *, '(a)' ) &
              '  Enter the spatial dimension M with the command'
            write ( *, '(a)' ) '    "M = ???"'
            exit
          end if

          call radius_make_2 ( m, a, b, n, cvt_m, cvt_center, &
            cvt_radius )
          exit
        end if

      end do

      call cvt_make ( alpha, beta, m, a, b, maxit, ns, density_function, &
        n, cvt_center )

    else if ( s_eqi ( command, 'CVT_OVERLAP' ) ) then

      call basis_overlap ( m, n, cvt_center, cvt_radius )

    else if ( s_eqi ( command, 'CVT_PLOT' ) ) then

      call basis_plot ( m, n, cvt_center, cvt_radius )

    else if ( s_eqi ( command, 'CVT_READ' ) ) then

      call cvt_read ( m, n, cvt_center, cvt_radius, cvt_file )

    else if ( s_eqi ( command, 'CVT_SCALE' ) ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Enter a scale factor to apply to the radius'
      write ( *, '(a)' ) 'of the current CVT function set:'

      read ( *, *, iostat = ios ) cvt_radius_factor
      if ( ios /= 0 ) then
        exit
      end if

      cvt_radius(1:n) = cvt_radius_factor * cvt_radius(1:n)

    else if ( s_eqi ( command, 'CVT_WRITE' ) ) then

      if ( n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Use the command cvt_make first!'
        cycle
      end if

      call cvt_write ( m, n, cvt_center, cvt_radius, cvt_file )

    else if ( s_eqi ( command(1:17), 'DENSITY_FUNCTION=' ) ) then

      call s_to_i4 ( command(18:), i_temp, ierror, last )

      if ( ierror == 0 ) then

        density_function = i_temp

        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  DENSITY_FUNCTION has been set to ', &
          density_function

      end if

    else if ( s_eqi ( command, 'HALTON_MAKE' ) ) then

      if ( m <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'HALTON_MAKE - Error!'
        write ( *, '(a)' ) '  The spatial dimension must be positive!'
        write ( *, '(a)' ) '  Enter the spatial dimension M with the command'
        write ( *, '(a)' ) '    "M = ???"'
        cycle
      end if

      if ( n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'HALTON_MAKE - Error!'
        write ( *, '(a)' ) '  The number of points N must be positive!'
        write ( *, '(a)' ) '  Enter the number of points with the command'
        write ( *, '(a)' ) '    "N = ???"'
        cycle
      end if

      call halton_generate ( m, n, skip, base, halton_points )

    else if ( s_eqi ( command, 'HALTON_READ' ) ) then

      call halton_read ( m, n, halton_points, halton_file )

    else if ( s_eqi ( command, 'HALTON_WRITE' ) ) then

      if ( n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'MESHLESS - Error!'
        write ( *, '(a)' ) '  No Halton basis functions have been defined.'
        write ( *, '(a)' ) '  Use the command HALTON_MAKE or HALTON_READ first!'
        cycle
      end if

      call halton_write ( m, n, skip, base, halton_points, halton_file )

    else if ( s_eqi ( command, 'HELP' ) ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Commands:'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  DENSITY_FUNCTION = ... specify density.'
      write ( *, '(a)' ) '  MAXIT = ...   specify number of iterations.'
      write ( *, '(a)' ) '  N =     ...   specify the number of basis points.'
      write ( *, '(a)' ) '  M =     ...   specify the spatial dimension.'
      write ( *, '(a)' ) '  NS =    ...   specify the sampling point density'
      write ( *, '(a)' ) '                for estimating areas.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  BASIS_MAKE    create a set of basis'
      write ( *, '(a)' ) '                functions.'
      write ( *, '(a)' ) '  BASIS_OVERLAP count the basis function overlaps.'
      write ( *, '(a)' ) '  BASIS_PLOT    plot a set of basis functions.'
      write ( *, '(a)' ) '  BASIS_READ    read a set of basis'
      write ( *, '(a)' ) '                functions from a file.'
      write ( *, '(a)' ) '  BASIS_SCALE   make the basis functions "wider"'
      write ( *, '(a)' ) '  BASIS_WRITE   write the current basis functions '
      write ( *, '(a)' ) '                to a file.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  CENTER_PLOT'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  CVT_MAKE      create a set of CVT basis'
      write ( *, '(a)' ) '                functions.'
      write ( *, '(a)' ) '  CVT_OVERLAP   count the basis function overlaps.'
      write ( *, '(a)' ) '  CVT_PLOT      plot a set of CVT basis functions.'
      write ( *, '(a)' ) '  CVT_READ      read a set of CVT basis'
      write ( *, '(a)' ) '                functions from a file.'
      write ( *, '(a)' ) '  CVT_SCALE     make CVT functions "wider")'
      write ( *, '(a)' ) '  CVT_WRITE     write the current CVT basis'
      write ( *, '(a)' ) '                functions to a file.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  HALTON_MAKE   create a set of Halton basis'
      write ( *, '(a)' ) '                functions.'
      write ( *, '(a)' ) '                Specify M and N first!'
      write ( *, '(a)' ) '  HALTON_READ   read a set of Halton basis'
      write ( *, '(a)' ) '                functions from a file.'
      write ( *, '(a)' ) '  HALTON_WRITE  write the current Halton basis '
      write ( *, '(a)' ) '                functions to a file.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  UNIFORM_MAKE  create a set of uniform basis'
      write ( *, '(a)' ) '                functions.'
      write ( *, '(a)' ) '  UNIFORM_READ  read a set of uniform basis'
      write ( *, '(a)' ) '                functions from a file.'
      write ( *, '(a)' ) '  UNIFORM_WRITE write the current uniform basis'
      write ( *, '(a)' ) '                functions to a file.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Help          print the list of commands.'
      write ( *, '(a)' ) '  Print         print the value of some quantity.'
      write ( *, '(a)' ) '  Quit          terminate the program.'

    else if ( s_eqi ( command(1:2), 'M=' ) ) then

      call s_to_i4 ( command(3:), i_temp, ierror, last )

      if ( ierror /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'MESHLESS - Error!'
        write ( *, '(a)' ) '  Could not set M!'
        cycle
      end if

      m = i_temp
      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) '  M has been set to ', m

      if ( allocated ( a ) ) then
        deallocate ( a ) 
      end if
      allocate ( a(m) )
      a(1:m) = -1.0D+00

      if ( allocated ( b ) ) then
        deallocate ( b ) 
      end if
      allocate ( b(m) )
      b(1:m) = +1.0D+00

      if ( allocated ( base ) ) then
        deallocate ( base )
      end if
      allocate ( base(m) )

      if ( allocated ( basis_center ) ) then
        deallocate ( basis_center )
      end if
      if ( 0 < n ) then
        allocate ( basis_center(m,n) )
      end if

      if ( allocated ( cvt_center ) ) then
        deallocate ( cvt_center )
      end if
      if ( 0 < n ) then
        allocate ( cvt_center(m,n) )
      end if

      if ( allocated ( halton_points ) ) then
        deallocate ( halton_points )
      end if
      if ( 0 < n ) then
        allocate ( halton_points(m,n) )
      end if

      if ( allocated ( uniform_points ) ) then
        deallocate ( uniform_points )
      end if
      if ( 0 < n ) then
        allocate ( uniform_points(m,n) )
      end if