meshless.f90

meshles是不采用网格

!
!  We do not allow the case where a point in CENTER is exactly equal to PT.
!
    if ( all ( pt(1:m) == center(1:m,i) ) ) then
      cycle
    end if
!
!  A center point is a suitable neighbor if its L-Infinity distance
!  from PT is no more than R.  But for some reason, we then record
!  its L2 distance...Is this an inconsistency?
!
    if ( all ( abs ( pt(1:m) - center(1:m,i) ) <= r ) ) then
      np = np + 1
      dist(np) = sqrt ( sum ( ( pt(1:m) - center(1:m,i) )**2 ) )
      ord(np) = i
    end if

  end do

  if ( np == 0 ) then
    ierr = -1
  else
    ierr = 0
  end if

  return
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
  implicit none

  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine halton_generate ( m, n, skip, base, r )

!*****************************************************************************80
!
!! HALTON_GENERATE generates a Halton dataset.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    16 March 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points to generate.
!
!    Input, integer SKIP, the number of initial points to skip.
!
!    Output, integer BASE(M), the bases used for the sequence.
!
!    Output, real ( kind = 8 ) R(M,N), the points.
!
  implicit none

  integer m
  integer n

  integer, dimension ( m ) :: base
  integer i
  integer j
  integer prime
  real ( kind = 8 ), dimension ( m, n ) :: r
  integer seed
  integer skip

  do i = 1, m
    base(i) = prime(i)
  end do

  do j = 1, n
    seed = skip + j - 1
    call i4_to_halton ( seed, base, m, r(1:m,j) )
  end do

  return
end
subroutine halton_read ( m, n, points, input_file )

!*****************************************************************************80
!
!! HALTON_READ reads Halton data from a file.
!
!  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 points.
!
!    Output, real ( kind = 8 ) POINTS(M,N), the Halton points.
!
!    Input, character ( len = * ) INPUT_FILE, the name of the file.
!
  implicit none

  integer n
  integer m

  integer i
  character ( len = * ) input_file
  real ( kind = 8 ), dimension ( m, n ) :: points

  write ( *, '(a)' ) ' '
  write ( *, '(a)' )  'HALTON_READ:'
  write ( *, '(a)' ) '  Reading Halton data file: ' // trim ( input_file )
  write ( *, '(a,i6)' ) '  Number of data points is ', n

  open ( unit = 12, file = input_file, form = 'formatted', status = 'old' )

  do i = 1, n
    read ( 12, * ) points(1:m,i)
  end do

  close ( unit = 12 )

  return
end
subroutine halton_write ( m, n, skip, base, r, file_out_name )

!*****************************************************************************80
!
!! HALTON_WRITE writes a Halton dataset to a file.
!
!  Discussion:
!
!    The initial lines of the file are comments, which begin with a
!    '#' character.
!
!    Thereafter, each line of the file contains the M-dimensional
!    components of the SKIP+I-1 entry of the Halton sequence.
!
!    For the Halton sequence, the value of SKIP is the same
!    as the value of SEED used to generate the first point.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    04 October 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of (successive) points.
!
!    Input, integer SKIP, the number of skipped points.
!
!    Input, integer BASE(M), the bases used for the sequence.
!
!    Input, real ( kind = 8 ) R(M,N), the points.
!
!    Input, character ( len = * ) FILE_OUT_NAME, the name of
!    the output file.
!
  implicit none

  integer m
  integer n

  integer base(m)
  character ( len = * ) file_out_name
  integer file_out_unit
  integer ios
  integer j
  integer mhi
  integer mlo
  real ( kind = 8 ) r(m,n)
  integer skip
  character ( len = 40 ) string

  call get_unit ( file_out_unit )

  open ( unit = file_out_unit, file = file_out_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HALTON_WRITE - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output file.'
    stop
  end if

  call timestring ( string )

  write ( file_out_unit, '(a)'      ) '#  ' // trim ( file_out_name )
  write ( file_out_unit, '(a)'      ) '#  created by MESHLESS.F90'
  write ( file_out_unit, '(a)'      ) '#'
  write ( file_out_unit, '(a)'      ) '#  File generated on ' // trim ( string )
  write ( file_out_unit, '(a)'      ) '#'
  write ( file_out_unit, '(a,i6)'   ) '#  Spatial dimension M = ', m
  write ( file_out_unit, '(a,i6)'   ) '#  Number of points N = ', n
  do mlo = 1, m, 10
    mhi = min ( mlo + 9, m )
    if ( mlo == 1 ) then
      write ( file_out_unit, '(a,20i6)' ) '#  Bases: ', base(mlo:mhi)
    else
      write ( file_out_unit, '(a,20i6)' ) '#         ', base(mlo:mhi)
    end if
  end do
  write ( file_out_unit, '(a,i6)'   ) '#  Initial values skipped = ', skip
  write ( file_out_unit, '(a)'      ) '#'

  write ( string, '(a,i3,a)' ) '(', m, 'f10.6)'

  do j = 1, n
    write ( file_out_unit, string ) r(1:m,j)
  end do

  close ( unit = file_out_unit )

  return
end
subroutine i4_swap ( i, j )

!*****************************************************************************80
!
!! I4_SWAP swaps two integer values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On output, the values of I and
!    J have been interchanged.
!
  implicit none

  integer i
  integer j
  integer k

  k = i
  i = j
  j = k

  return
end
subroutine i4_to_halton ( seed, base, ndim, r )

!*****************************************************************************80
!
!! I4_TO_HALTON computes an element of a Halton sequence.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    26 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    John Halton,
!    On the efficiency of certain quasi-random sequences of points
!    in evaluating multi-dimensional integrals,
!    Numerische Mathematik,
!    Volume 2, pages 84-90, 1960.
!
!  Parameters:
!
!    Input, integer SEED, the index of the desired element.
!    Only the absolute value of SEED is considered.  SEED = 0 is allowed,
!    and returns R = 0.
!
!    Input, integer BASE(NDIM), the Halton bases, which should be
!    distinct prime numbers.  This routine only checks that each base
!    is greater than 1.
!
!    Input, integer NDIM, the dimension of the sequence.
!
!    Output, real ( kind = 8 ) R(NDIM), the SEED-th element of the
!    Halton sequence for the given bases.
!
  implicit none

  integer ndim

  integer base(ndim)
  real ( kind = 8 ) base_inv(ndim)
  integer digit(ndim)
  integer i
  real ( kind = 8 ) r(ndim)
  integer seed
  integer seed2(ndim)

  seed2(1:ndim) = abs ( seed )

  r(1:ndim) = 0.0D+00

  if ( any ( base(1:ndim) <= 1 ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I4_TO_HALTON - Fatal error!'
    write ( *, '(a)' ) '  An input base BASE is <= 1!'
    do i = 1, ndim
      write ( *, '(i6,i6)' ) i, base(i)
    end do
    stop
  end if

  base_inv(1:ndim) = 1.0D+00 / real ( base(1:ndim), kind = 8 )

  do while ( any ( seed2(1:ndim) /= 0 ) )
    digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim) )
    r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), kind = 8 ) * base_inv(1:ndim)
    base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), kind = 8 )
    seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
  end do

  return
end
function point_inside_box_2d ( x1, y1, x2, y2, x, y )

!*****************************************************************************80
!
!! POINT_INSIDE_BOX_2D determines if a point is inside a box in 2D.
!
!  Discussion:
!
!    A "box" is defined by its "left down" corner and its
!    "right up" corner, and all the points between.  It is
!    assumed that the sides of the box align with coordinate directions.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X1, Y1, X2, Y2, the two corners of the box.
!
!    Input, real ( kind = 8 ) X, Y, the point to be checked.
!
!    Output, logical POINT_INSIDE_BOX_2D, is .TRUE. if (X,Y) is inside the
!    box, or on its boundary, and .FALSE. otherwise.
!
  implicit none

  logical point_inside_box_2d
  real ( kind = 8 ) x
  real ( kind = 8 ) x1
  real ( kind = 8 ) x2
  real ( kind = 8 ) y
  real ( kind = 8 ) y1
  real ( kind = 8 ) y2

  if ( x1 <= x .and. x <= x2 .and. &
       y1 <= y .and. y <= y2 ) then
    point_inside_box_2d = .true.
  else
    point_inside_box_2d = .false.
  end if

  return
end
function prime ( n )

!*****************************************************************************80
!
!! PRIME returns any of the first PRIME_MAX prime numbers.
!
!  Discussion:
!
!    PRIME_MAX is 1500, and the largest prime stored is 12553.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    21 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964, pages 870-873.
!
!    Daniel Zwillinger,
!    CRC Standard Mathematical Tables and Formulae,
!    30th Edition,
!    CRC Press, 1996, pages 95-98.
!
!  Parameters:
!
!    Input, integer N, the index of the desired prime number.
!    N = -1 returns PRIME_MAX, the index of the largest prime available.
!    N = 0 is legal, returning PRIME = 1.
!    It should generally be true that 0 <= N <= PRIME_MAX.
!
!    Output, integer PRIME, the N-th prime.  If N is out of range, PRIME
!    is returned as 0.
!
  implicit none

  integer, parameter :: prime_max = 1500

  integer, save :: icall = 0
  integer n
  integer, save, dimension ( prime_max ) :: npvec
  integer prime

  if ( icall == 0 ) then

    icall = 1

    npvec(1:100) = (/ &
        2,    3,    5,    7,   11,   13,   17,   19,   23,   29, &
       31,   37,   41,   43,   47,   53,   59,   61,   67,   71, &
       73,   79,   83,   89,   97,  101,  103,  107,  109,  113, &
      127,  131,  137,  139,  149,  151,  157,  163,  167,  173, &
      179,  181,  191,  193,  197,  199,  211,  223,  227,  229, &
      233,  239,  241,  251,  257,  263,  269,  271,  277,  281, &
      283,  293,  307,  311,  313,  317,  331,  337,  347,  349, &
      353,  359,  367,  373,  379,  383,  389,  397,  401,  409, &
      419,  421,  431,  433,  439,  443,  449,  457,  461,  463, &
      467,  479,  487,  491,  499,  503,  509,  521,  523,  541 /)

    npvec(101:200) = (/ &
      547,  557,  563,  569,  571,  577,  587,  593,  599,  601, &
      607,  613,  617,  619,  631,  641,  64