meshless.f90

meshles是不采用网格

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

  n = 0

  do

    read ( 12, *, iostat = ios ) c_temp(1:m), r_temp

    if ( ios /= 0 ) then
      exit
    end if

    n = n + 1
    center(1:m,n) = c_temp(1:m)
    radius(n) = r_temp

  end do

  close ( unit = 12 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BASIS_READ:'
  write ( *, '(a,i6)' ) '  Number of data points read was ', n

  return
end
subroutine basis_write ( m, n, center, radius, output_file )

!*****************************************************************************80
!
!! BASIS_WRITE writes the basis data to 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.
!
!    Input, real ( kind = 8 ) CENTER(M,N), the basis function centers.
!
!    Input, real ( kind = 8 ) RADIUS(N), the basis function radii.
!
!    Input, character ( len = * ) OUTPUT_FILE, the name of the output file.
!
  implicit none

  integer n
  integer m

  integer i
  character ( len = * ) output_file
  real ( kind = 8 ), dimension ( m, n ) :: center
  real ( kind = 8 ), dimension ( n ) :: radius

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BASIS_WRITE:'
  write ( *, '(a)' ) '  Write basis data file: ' // trim ( output_file )
  write ( *, '(a,i6)' ) '  Number of data points is ', n

  open ( unit = 12, file = output_file, form = 'formatted', status = 'replace' )

  do i = 1, n
    write ( 12, * ) center(1:m,i), radius(i)
  end do

  close ( unit = 12 )

  return
end
subroutine ch_cap ( c )

!*****************************************************************************80
!
!! CH_CAP capitalizes a single character.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none

  character c
  integer itemp

  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
subroutine center_plot ( m, n, center, a, b )

!*****************************************************************************80
!
!! CENTER_PLOT plots the center in the region.
!
!  Discussion:
!
!    My idea of plotting the basis functions is OK if there are just a
!    few, but quickly becomes useless.  Here, at least, I can plot just
!    the center points, and get an idea of the distribution.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    24 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 ) 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 ) blue
  real ( kind = 8 ) center(m,n)
  character ( len = 80 ) file_name
  logical, parameter :: filled = .false.
  real ( kind = 8 ) green
  integer ierror
  integer iunit
  integer marker_size
  integer nx
  integer ny
  real ( kind = 8 ) red
  real ( kind = 8 ) rxmax
  real ( kind = 8 ) rxmin
  real ( kind = 8 ) rymax
  real ( kind = 8 ) rymin
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
!
!  NOT CORRECT FOR 2 < M.
!
  rxmin = a(1)
  rxmax = b(1)
  rymin = a(2)
  rymax = b(2)

  iunit = 1
  file_name = 'center_plot.ps'

  call ps_file_open ( file_name, iunit, ierror )

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

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

  call ps_page_head ( xmax, xmin, ymax, ymin )
!
!  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.
!
  red = 0.2D+00
  green = 0.3D+00
  blue = 0.4D+00

  call ps_color_line_set ( red, green, blue )

  nx = 11
  ny = 11

  call ps_grid_cartesian ( rxmin, rxmax, nx, rymin, rymax, ny )

  marker_size = 3
  call ps_marker_size ( marker_size )
!
!  Draw the center points.
!
  x(1:n) = center(1,1:n)
  y(1:n) = center(2,1:n)

  call ps_mark_circles ( n, x, y )

  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

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

  return
end
subroutine cvt_make ( alpha, beta, m, a, b, maxit, ns, density_function, &
  n, center )

!*****************************************************************************80
!
!! CVT_MAKE computes centroidal Voronoi tessellation generators.
!
!  Discussion:
!
!    The routine is given an initial set of points that define a 
!    Voronoi tessellation.  It computes new points that define a Voronoi
!    tessellation, and which have the property that they are the centroids
!    of their Voronoi regions.
!
!    No stopping criterion is used, for now, except to iterate to MAXIT.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    29 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ALPHA(2), BETA(2), coefficients used in the
!    generation program.  ALPHA(1) and ALPHA(2) must lie between 0 and 1
!    and sum to 1.  The same is true for BETA.  Common values are
!    ALPHA(1) = BETA(1) = 0,
!    ALPHA(2) = BETA(2) = 1.
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = 8 ) A(M), B(M), the coordinates of the
!    two extreme corners of the box that defines the region.
!
!    Input, integer MAXIT, the maximum number of iterations.
!
!    Input, integer NS, the average number of sampling points 
!    per generator, on each step of the iteration.
!
!    Input, integer DENSITY_FUNCTION, specifies 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) ) )
!
!    Input, integer N, the number of generators.
!
!    Input/output, real ( kind = 8 ) CENTER(M,N).
!    On input, initial points that generate the Voronoi regions.
!    On output, points that generate the Voronoi regions, which are
!    also the centroids of those regions.
!
!  Local variables:
!
!    Integer UPDATES(N), counts the number of times a generator 
!    has been updated.
!
  implicit none

  integer n
  integer m

  real ( kind = 8 ) a(1:m)
  real ( kind = 8 ) alpha(2)
  real ( kind = 8 ) b(1:m)
  real ( kind = 8 ) beta(2)
  real ( kind = 8 ) center(m,n)
  real ( kind = 8 ) center2(m,n)
  integer count(n)
  integer density_function
  integer ic
  integer iloop
  integer j
  integer k
  integer maxit
  integer ns
  real ( kind = 8 ) p1(m)
  real ( kind = 8 ) p2(m)
  real ( kind = 8 ) s
  integer updates(n)
  real ( kind = 8 ) y(m)

  if ( m <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CVT_MAKE - Error!'
    write ( *, '(a)' ) '  The spatial dimension must be positive!'
    write ( *, '(a)' ) '  Enter a spatial dimension with the "M = " command.'
    return
  end if
!
!  UPDATES starts at 1.
!
  updates(1:n) = 1
!
!  Iterate MAXIT times.
!
  do iloop = 1, maxit

    center2(1:m,1:n) = 0.0D+00
    count(1:n) = 0

    do j = 1, n * ns
!
!  Generate a random sampling point Y.
!
      call random_generator ( m, a, b, density_function, y )
!
!  Find the nearest generator CENTER.  Its index is IC.
!
      call find_closest ( m, y, n, center, ic, s )
!
!  Add Y to the averaging data for CENTER(*,IC).
!
      center2(1:m,ic) = center2(1:m,ic) + y(1:m)
      count(ic) = count(ic) + 1

    end do
!
!  For each cell J, average the estimated centroid CENTER with the 
!  averaged hit points CENTER2.  (There are COUNT(J) of these new points).
!
    do j = 1, n

      if ( count(j) /= 0.0 ) then
!
!  This loop and the next could be combined, eliminating the
!  temporary variables P1 and P2.
!
        do k = 1, m
          
          p1(k) = ( alpha(1) * updates(j) + beta(1) ) * center(k,j)
          p2(k) = ( alpha(2) * updates(j) + beta(2) ) * center2(k,j) &
            / dble ( count(j) )
        end do

        do k = 1, m
          center(k,j) = ( p1(k) + p2(k) ) / dble ( updates(j) + 1 )
        end do

        updates(j) = updates(j) + 1

      end if

    end do

  end do

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

!*****************************************************************************80
!
!! CVT_READ reads the CVT 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 ) CENTER(M,N), the centroidal 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 ) :: center
  real ( kind = 8 ), dimension ( n ) :: radius

  write ( *, '(a)' ) ' '
  write ( *, '(a)')  'CVT_READ:'
  write ( *, '(a)' ) '  Reading Centroidal Voronoi 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, * ) center(1:m,i), radius(i)
  end do

  close ( unit = 12 )

  return
end
subroutine cvt_write ( m, n, center, radius, output_file )

!*****************************************************************************80
!
!! CVT_WRITE writes the CVT data to a file.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    05 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 centroidal points.
!
!    Input, character ( len = * ) OUTPUT_FILE, the name of the output file.
!
  implicit none

  integer n
  integer m

  integer i
  character ( len = * ) output_file
  real ( kind = 8 ), dimension ( m, n ) :: center
  real ( kind = 8 ), dimension ( n ) :: radius

  write ( *, '(a)' ) ' '
  write ( *, '(a)' )  'CVT_WRITE:'
  write ( *, '(a)' ) '  Write Centroidal Voronoi data file: ' &
    // trim ( output_file )
  write ( *, '(a,i6)' ) '  Number of data points is ', n

  open ( unit = 12, file = output_file, form = 'formatted', status = 'replace' )

  do i = 1, n
    write ( 12, * ) center(1:m,i), radius(i)
  end do

  close ( unit = 12 )

  return
end
function density ( m, x, density_function )

!*****************************************************************************80
!
!! DENSITY evaluates the density function.
!
!  Discussion:
!
!    The density function is used to control the generation of random points
!    that sample the region.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.