meshless.f90

meshles是不采用网格

!
!  Modified:
!
!    19 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = 8 ) X(M), the location of the point.
!
!    Input, integer DENSITY_FUNCTION, chooses the density function.
!    1: d(x) = 1.0;
!    2: d(x) = exp ( - 4.0 * ( sum(x(1:n)**2) ) )
!    3: d(x) = exp ( - 3.0 * ( 1.0 - sum(x(1:n)**2) ) )
!
!    Output, real ( kind = 8 ) DENSITY, the value of the density 
!    function, which should be between 0 and 1 in the region.
!
  implicit none

  integer m

  real ( kind = 8 ) density
  integer density_function
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) x(m)
!
!  Here is an odd feature that sets density to zero if the point
!  is more than 1 away from the origin.  This is not what I documented,
!  and should probably be controllable by the user!
!
  if ( 1.0D+00 < sum ( x(1:m)**2 ) ) then
    density = 0.0D+00
    return
  end if

  if ( density_function == 1 ) then
    density = 1.0D+00
  else if ( density_function == 2 ) then
    density = exp ( - 4.0D+00 * sum ( x(1:m)**2 ) )
  else if ( density_function == 3 ) then
    density = exp ( - 3.0D+00 * ( 1.0D+00 - sum ( x(1:m)**2 ) ) )
  else
    density = 1.0D+00
  end if

  return
end
subroutine eps_file_head ( file_name, x_ps_min, y_ps_min, x_ps_max, &
  y_ps_max )

!*****************************************************************************80
!
!! EPS_FILE_HEAD writes header information to an encapsulated PostScript file.
!
!  Discussion:
!
!    The file should contain the description of only one page, but this
!    is not currently checked.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 April 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Henry McGilton and Mary Campione,
!    PostScript by Example,
!    Addison-Wesley,
!    ISBN: 0-201-63228-4
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the output file.
!
!    Input, integer X_PS_MIN, Y_PS_MIN, X_PS_MAX, Y_PS_MAX, the minimum 
!    and maximum X and Y values of the data, in PostScript units.  Any data
!    that lies outside this range will not show up properly.  A reasonable
!    set of values might be 0, 0, 612, 792, or, for a half inch margin,
!    36, 36, 576, 756.
!
  implicit none

  character ( len = 8 ) date
  character ( len = * ) file_name
  real ( kind = 8 ) line_blue
  real ( kind = 8 ) line_green
  real ( kind = 8 ) line_red
  integer state
  integer unit
  integer x_ps_max
  integer x_ps_min
  integer y_ps_max
  integer y_ps_min
!
!  Determine if the PostScript state is acceptable.
!
  call ps_setting_int ( 'GET', 'STATE', state )

  if ( state /= 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EPS_FILE_HEAD - Fatal error!'
    write ( *, '(a,i9)' ) '  PostScript state is ', state
    write ( *, '(a)' ) '  PostScript state 1 is required.'
    return
  end if
!
!  Initialization
!
  call ps_default ( )
!
!  Get the unit number.
!
  call ps_setting_int ( 'GET', 'UNIT', unit )

  call date_and_time ( date )
!
!  Write the prolog.
!
  write ( unit, '(a)' )     '%!PS-Adobe-3.0 EPSF-3.0'
  write ( unit, '(a)' )     '%%Creator: ps_write.f90'
  write ( unit, '(a)' )     '%%Title: ' // trim ( file_name )
  write ( unit, '(a)' )     '%%CreationDate: '// trim ( date )
  write ( unit, '(a)' )     '%%Pages: 1'
  write ( unit, '(a,4i6)' ) '%%BoundingBox:', &
    x_ps_min, y_ps_min, x_ps_max, y_ps_max
  write ( unit, '(a)' )     '%%Document-Fonts: Times-Roman'
  write ( unit, '(a)' )     '%%LanguageLevel: 1'
  write ( unit, '(a)' )     '%%EndComments'
  write ( unit, '(a)' )     '%%BeginProlog'
  write ( unit, '(a)' )     '/inch {72 mul} def'
  write ( unit, '(a)' )     '%%EndProlog'
!
!  Set the font.
!
  write ( unit, '(a)' ) '/Times-Roman findfont'
  write ( unit, '(a)' ) '1.00 inch scalefont'
  write ( unit, '(a)' ) 'setfont'
!
!  Set the line color.
!
  line_red = 0.0D+00
  line_green = 0.0D+00
  line_blue = 0.0D+00

  call ps_color_line ( 'SET', line_red, line_green, line_blue )
!
!  Reset the state.
!
  state = 2

  call ps_setting_int ( 'SET', 'STATE', state )

  return
end
subroutine eps_file_tail ( )

!*****************************************************************************80
!
!! EPS_FILE_TAIL writes trailer information to an encapsulated PostScript file.
!
!  Discussion:
!
!    Looks like that penultimate 'end' line is not wanted, so I commented
!    it out.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    05 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Henry McGilton and Mary Campione,
!    PostScript by Example,
!    Addison-Wesley,
!    ISBN: 0-201-63228-4
!
!  Parameters:
!
!    None
!
  implicit none

  integer num_pages
  integer state
  integer unit
!
!  Determine if the PostScript state is acceptable.
!
  call ps_setting_int ( 'GET', 'STATE', state )

  if ( state == 3 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EPS_FILE_TAIL - Warning!'
    write ( *, '(a)' ) '  A page was open.  It is being forced closed.'
    state = 2
    call ps_setting_int ( 'SET', 'STATE', state )
  end if

  if ( state /= 2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EPS_FILE_TAIL - Fatal error!'
    write ( *, '(a,i9)' ) '  PostScript state is ', state
    write ( *, '(a)' ) '  PostScript state 2 is required.'
    return
  end if
!
!  Get the unit number.
!
  call ps_setting_int ( 'GET', 'UNIT', unit )
!
!  Retrieve the number of pages.
!
  call ps_setting_int ( 'GET', 'NUM_PAGES', num_pages )

  if ( 1 < num_pages ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EPS_FILE_TAIL - Warning!'
    write ( *, '(a)' ) '  An encapsulated PostScript file describes ONE page.'
    write ( *, '(a,i9,a)' ) '  This file describes ', num_pages, ' pages.'
    write ( *, '(a)' ) '  It is not a legal EPS file.'
  end if
!
!  Write the epilog.
!
  write ( unit, '(a)' ) '%%Trailer'
! write ( unit, '(a)' ) 'end'
  write ( unit, '(a)' ) '%%EOF'
!
!  Zero out the number of pages.
!
  num_pages = 0

  call ps_setting_int ( 'SET', 'NUM_PAGES', num_pages )
!
!  Reset the state.
!
  state = 4

  call ps_setting_int ( 'SET', 'STATE', state )

  return
end
subroutine file_column_count ( file_name, ncolumn )

!*****************************************************************************80
!
!! FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!    Each line of the file is presumed to consist of NCOLUMN words, separated
!    by spaces.
!
!    The routine reads the first line and counts the number of words.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NCOLUMN, the number of columns assumed to be in the file.
!
  implicit none

  character ( len = * ) file_name
  integer ios
  integer iunit
  character ( len = 256 ) line
  integer ncolumn
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    ncolumn = - 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file:'
    write ( *, '(a)' ) '    ' // trim ( file_name )
    return
  end if
!
!  Read one line.
!
  read ( iunit, '(a)', iostat = ios ) line

  if ( ios /= 0 ) then
    ncolumn = -2
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, '(a)' ) '  The initial line of the file could not be read.'
    write ( *, '(a)' ) '    ' // trim ( file_name )
    return
  end if

  close ( unit = iunit )

  call s_word_count ( line, ncolumn )

  return
end
subroutine file_line_count ( file_name, nline )

!*****************************************************************************80
!
!! FILE_LINE_COUNT counts the number of lines in a file.
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    21 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NLINE, the number of lines found in the file.
!
  implicit none

  character ( len = * ) file_name
  integer ios
  integer iunit
  integer nline

  nline = 0
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    nline = - 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file:'
    write ( *, '(a)' ) '    ' // trim ( file_name )
    return
  end if
!
!  Count the lines.
!
  do

    read ( iunit, '(a)', iostat = ios )

    if ( ios /= 0 ) then
      exit
    end if

    nline = nline + 1

  end do

  close ( unit = iunit )

  return
end
subroutine find_closest ( m, y, n, center, ic, s )

!*****************************************************************************80
!
!! FIND_CLOSEST finds the center point CENTER(1:M,IC) closest to Y(1:M).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = 8 ) Y(M), the point to be checked.
!
!    Input, integer N, the number of center points.
!
!    Input, real ( kind = 8 ) CENTER(M,N), the center points.
!
!    Output, integer IC, the index of the nearest center point.
!
!    Output, real ( kind = 8 ) S, the distance.
!
  implicit none

  integer n
  integer m

  real ( kind = 8 ) center(m,n)
  real ( kind = 8 ) dist_sq
  integer i
  integer ic
  real ( kind = 8 ) s
  real ( kind = 8 ) y(m)

  s = huge ( s )

  do i = 1, n

    dist_sq = sum ( ( center(1:m,i) - y(1:m) )**2 )

    if ( dist_sq < s ) then
      s = dist_sq
      ic = i
    end if

  end do

  s = sqrt ( s )

  return
end
subroutine find_re ( m, pt, r, n, center, np, dist, ord, ierr )

!*****************************************************************************80
!
!! FIND_RE finds the number of center points near a given point.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = 8 ) PT(M), the base point.
!
!    Input, real ( kind = 8 ) R, the maximum distance from PT 
!    to be considered.
!
!    Input, integer N, the number of center points.
!
!    Input, real ( kind = 8 ) CENTER(M,N), the center points.
!
!    Output, integer NP, the number of neighboring center points found.
!
!    Output, real ( kind = 8 ) DIST(1:NP), lists the Euclidean distance 
!    between PT and each of the neighboring center points found.
!
!    Output, integer ORD(1:NP), lists the indices in CENTER of the points found.
!
!    Output, integer IERR, error flag.
!    -1, if no points at all were found.
!    0, if at least one point was found.
!
  implicit none

  integer n
  integer m

  real ( kind = 8 ) center(m,n)
  real ( kind = 8 ) dist(n)
  integer i
  integer ierr
  integer np
  integer ord(n)
  real ( kind = 8 ) pt(m)
  real ( kind = 8 ) r
!
!  Search for a center point CENTER(I) which is no more than R away in X or Y,
!  and which is not identical to PT.  Count the number of such points found,
!  their distance, and indexes.
!
  np = 0

  do i = 1, n