meshless.f90

meshles是不采用网格

    else if ( s_eqi ( command(1:6), 'MAXIT=' ) ) then

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

      if ( ierror == 0 ) then

        maxit = i_temp
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  MAXIT has been set to ', maxit

      end if

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

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

      if ( ierror == 0 ) then

        n = i_temp
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  N has been set to ', n

      end if

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

      if ( allocated ( basis_radius ) ) then
        deallocate ( basis_radius )
      end if
      allocate ( basis_radius(n) )

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

      if ( allocated ( cvt_radius ) ) then
        deallocate ( cvt_radius )
      end if
      allocate ( cvt_radius(n) )

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

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

    else if ( s_eqi ( command(1:3), 'NS=' ) ) then

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

      if ( ierror == 0 ) then

        ns = i_temp
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  NS has been set to ', ns

      end if

    else if ( s_eqi ( command(1:5), 'PRINT' ) ) then

      what = adjustl ( command(6:) )
      if ( len_trim ( what ) == 0 ) then
        what = 'ALL'
      end if

      if ( s_eqi ( what, 'N' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, '(a,i6)' ) '  The number of points, N = ', n
      end if

      if ( s_eqi ( what, 'M' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, '(a,i6)' ) '  The spatial dimension, M = ', m
      end if

      if ( s_eqi ( what, 'MAXIT' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, '(a,i6)' ) &
          '  The maximum number of iterations, MAXIT = ', maxit
      end if

      if ( s_eqi ( what, 'NS' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, '(a,i6)' ) '  I don''t know what this is, NS = ', ns
      end if

      if ( s_eqi ( what, 'DENSITY_FUNCTION' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, '(a,i6)' ) '  The density function, DENSITY_FUNCTION = ', &
          density_function
      end if

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

      exit

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

      if ( m <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'UNIFORM_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)' ) 'UNIFORM_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 uniform_make ( m, n, a, b, density_function, &
        uniform_points )

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

      call uniform_read ( m, n, uniform_points, uniform_file )

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

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

      call uniform_write ( m, n, uniform_points, uniform_file )

    else

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'MESHLESS - Warning'
      write ( *, '(a)' ) '  Your command was not recognized.'

    end if

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'MESHLESS'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine basis_box ( m, n, center, radius, a, b )

!*****************************************************************************80
!
!! BASIS_BOX finds a box that contains all the basis functions.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    26 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) CENTER(M,N), the basis function centers.
!
!    Input, real ( kind = 8 ) RADIUS(N), the basis function radii.
!
!    Output, real ( kind = 8 ) A(M), B(M), the minimum and maximum
!    values in each dimension.
!
  implicit none

  integer n
  integer m

  real ( kind = 8 ) a(m)
  real ( kind = 8 ) b(m)
  real ( kind = 8 ), dimension ( m, n ) :: center
  integer i
  real ( kind = 8 ) radius(n)

  do i = 1, m
    a(i) = minval ( center(i,1:n) - radius(1:n) )
    b(i) = maxval ( center(i,1:n) + radius(1:n) )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BASIS_BOX'
  write ( *, '(a)' ) '  Compute the bounding box.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Box limits for center points:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Dimension  Lower  Upper'
  write ( *, '(a)' ) ' '

  do i = 1, m
    write ( *, '(i4,2g14.6)' ) &
      i, minval ( center(i,1:n) ), maxval ( center(i,1:n) )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Box limits for radius:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Lower  Upper'
  write ( *, '(a)' ) ' '
  write ( *, '(2g14.6)' ) minval ( radius(1:n) ), maxval ( radius(1:n) )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Box limits for basis support:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Dimension  Lower  Upper'
  write ( *, '(a)' ) ' '

  do i = 1, m
    write ( *, '(i4,2g14.6)' ) i, a(i), b(i)
  end do

  return
end
subroutine basis_overlap ( m, n, center, radius )

!*****************************************************************************80
!
!! BASIS_OVERLAP counts the number of basis functions with overlapping support.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of basis functions.
!
!    Input, real ( kind = 8 ) CENTER(M,N), the basis function centers.
!
!    Input, real ( kind = 8 ) RADIUS(N), the basis function radii.
!
  implicit none

  integer n
  integer m

  real ( kind = 8 ), dimension ( m, n ) :: center
  real ( kind = 8 ) dist
  integer i
  integer j
  integer mel
  integer nel
  real ( kind = 8 ), dimension ( n ) :: radius

  if ( n <= 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BASIS_OVERLAP - Warning!'
    write ( *, '(a)' ) '  N <= 1, so no overlap possible.'
    return
  end if

  mel = ( n * ( n - 1 ) ) / 2
  nel = 0

  do i = 1, n
    do j = i+1, n

      dist = sqrt ( sum ( ( center(1:m,i) - center(1:m,j) )**2 ) )

      if ( dist <= ( radius(i) + radius(j) ) ) then
        nel = nel + 1
      end if

    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BASIS_OVERLAP:'
  write ( *, '(a)' ) '  Count the pairs of basis functions '
  write ( *, '(a)' ) '  with overlapping support.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Maximum possible count = ', mel
  write ( *, '(a,i6)' ) '  Actual count =           ', nel
  write ( *, '(a,f10.4)' ) &
    '  Percentage =             ', &
    real ( 100 * nel, kind = 8 ) / real ( mel, kind = 8 )

  return
end
subroutine basis_plot ( m, n, center, radius )

!*****************************************************************************80
!
!! BASIS_PLOT plots the basis functions in the region.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of basis functions.
!
!    Input, real ( kind = 8 ) CENTER(M,N), the basis function centers.
!
!    Input, real ( kind = 8 ) RADIUS(N), the basis function radii.
!
  implicit none

  integer n
  integer m

  real ( kind = 8 ) blue
  real ( kind = 8 ) center(m,n)
  character ( len = 80 ) file_name
  logical, parameter :: filled = .false.
  real ( kind = 8 ) green
  integer i
  integer ierror
  integer iunit
  integer nx
  integer ny
  real ( kind = 8 ) r
  real ( kind = 8 ) radius(n)
  real ( kind = 8 ) red
  real ( kind = 8 ) rxmax
  real ( kind = 8 ) rxmin
  real ( kind = 8 ) rymax
  real ( kind = 8 ) rymin
  real ( kind = 8 ) x
  integer, parameter :: x_ps_max = 576
  integer, parameter :: x_ps_min = 36
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) y
  integer, parameter :: y_ps_max = 756
  integer, parameter :: y_ps_min = 36
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin

  do i = 1, n
    if ( i == 1 ) then
      rxmin = center(1,i) - radius(i)
      rxmax = center(1,i) + radius(i)
      rymin = center(2,i) - radius(i)
      rymax = center(2,i) + radius(i)
    else
      rxmin = min ( rxmin, center(1,i) - radius(i) )
      rxmax = max ( rxmax, center(1,i) + radius(i) )
      rymin = min ( rymin, center(2,i) - radius(i) )
      rymax = max ( rymax, center(2,i) + radius(i) )
    end if
  end do

  write ( *, * ) 'RXMIN, RXMAX = ', rxmin, rxmax
  write ( *, * ) 'RYMIN, RYMAX = ', rymin, rymax

  call get_unit ( iunit )

  file_name = 'basis.eps'

  call ps_file_open ( file_name, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BASIS_PLOT'
    write ( *, '(a,i6)' ) '  File creation error ', ierror
    return
  end if

  call eps_file_head ( file_name, x_ps_min, y_ps_min, x_ps_max, &
    y_ps_max )

  xmin = rxmin - 0.1D+00 * ( rxmax - rxmin )
  xmax = rxmax + 0.1D+00 * ( rxmax - rxmin )
  ymin = rymin - 0.1D+00 * ( rymax - rymin )
  ymax = rymax + 0.1D+00 * ( rymax - rymin )

  call ps_page_head ( xmin, ymin, xmax, ymax )
!
!  Draw the outline of the region.
!
  red = 0.0D+00
  green = 0.0D+00
  blue = 1.0D+00

  call ps_color_line_set ( red, green, blue )

  call ps_moveto ( rxmin, rymin )
  call ps_lineto ( rxmax, rymin )
  call ps_lineto ( rxmax, rymax )
  call ps_lineto ( rxmin, rymax )
  call ps_lineto ( rxmin, rymin )
!
!  Draw a grid in the region.
!
  nx = 11
  ny = 11

  red = 0.5D+00
  green = 0.5D+00
  blue = 0.5D+00

  call ps_grid_cartesian ( rxmin, rxmax, nx, rymin, rymax, ny )
!
!  Draw the basis function support disks.
!
  red = 0.2D+00
  green = 0.2D+00
  blue = 0.2D+00

  call ps_color_line_set ( red, green, blue )

  red = 0.9D+00
  green = 0.6D+00
  blue = 0.6D+00

  call ps_color_fill_set ( red, green, blue )

  do i = 1, n

    x = center(1,i)
    y = center(2,i)
    r = radius(i)

    call ps_mark_disk ( x, y )

    if ( filled ) then
      call ps_circle_fill ( x, y, r )
    end if

    call ps_circle ( x, y, r )

  end do

  call ps_page_tail ( )

  call eps_file_tail ( )

  call ps_file_close ( iunit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BASIS_PLOT'
  write ( *, '(a)' ) &
    '  Created a basis function plot file ' // trim ( file_name )

  return
end
subroutine basis_read ( m, n, center, radius, input_file )

!*****************************************************************************80
!
!! BASIS_READ reads basis data from a file.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    24 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Output, integer N, the number of points.
!
!    Output, real ( kind = 8 ) CENTER(M,N), the basis function centers.
!
!    Output, real ( kind = 8 ) RADIUS(N), the basis function radii.
!
!    Input, character ( len = * ) INPUT_FILE, the name of the file.
!
  implicit none

  integer m

  real ( kind = 8 ), dimension ( m ) :: c_temp
  real ( kind = 8 ), dimension ( m, * ) :: center
  character ( len = * ) input_file
  integer ios
  integer n
  real ( kind = 8 ) r_temp
  real ( kind = 8 ), dimension ( * ) :: radius

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BASIS_READ:'
  write ( *, '(a)' ) '  Reading basis data file: ' // trim ( input_file )