2次元Delaunay法(杉原厚吉先生)

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88

 
 
 
-
|
|
|
|
!
 
 
 
 
 
 
-
|
|
!
 
 
 
 
-
|
|
!
 
 
 
 
-
|
|
!
 
 
 
 
 
 
-
!
 
 
-
!
 
 
 
-
!
 
 
 
 
 
-
!
 
 
-
!
 
 
-
!
 
 
-
!
-
!
 
 
 
 
-
|
|
!
 
 
 
-
!
 
 

      SUBROUTINE DRDELVTK(NOFP,IX,IY,IXYMAX,RADIUS,
     +  KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE,
     +  KVCOLO,KFCOLO)
* ------------------------------------------------------------ *
*                                                              *
*      Draw the Delaunay triangulation with VTK format         *
*                                                              *
* ------------------------------------------------------------ *
      DIMENSION IX(1),IY(1)
      DIMENSION KPARA(1),KVE(1),KFE(1),KSV(1),KEV(1),KRF(1),KLF(1)
      DIMENSION KSCE(1),KSCCE(1),KECE(1),KECCE(1)
      DIMENSION KVCOLO(1),KFCOLO(1)
      DIMENSION LISTV(100),LISTE(100)
      DIMENSION LSTUUF(KPARA(2)),LSTUF(KPARA(2))
      integer NOFVE
*
*  ---  open the file
*
      OPEN(80,file='Delaunay.vtk',FORM='FORMATTED',STATUS='UNKNOWN')
      write(80,"(a)") '# vtk DataFile Version 3.0'
      write(80,"(a)") 'vtk output'
      write(80,"(a)") 'ASCII'
      write(80,"(a)") 'DATASET UNSTRUCTURED_GRID'
*
*  ---  write the input coodinate
*
      write(80,"(a6,i8,a8)") 'POINTS', NOFP, 'float'
      do II=1,NOFP
        write(80,*) IX(II), IY(II), 0
      end do
 
*
*  ---  write connectivity
*
      NOFF=0
      KUF = KPARA(4)
      do while (KUF > 0)
        NOFF = NOFF + 1
        LSTUUF(NOFF) = KUF
        KUF = KFE(KUF)
      end do
*
      NOF=0
      do II=1,KPARA(2)
        if (KFCOLO(II) /= 1) cycle
*
        do JJ=1,NOFF
          if (LSTUUF(JJ) == II) exit
        end do        
        if (JJ <= NOFF) cycle
*
        call VEONF(II,
     +  KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE,
     +  NOFVE,LISTV,LISTE)
        AREA = (IX(LISTV(2))-IX(LISTV(1)))*(IY(LISTV(3))-IY(LISTV(1)))
     +       - (IX(LISTV(3))-IX(LISTV(1)))*(IY(LISTV(2))-IY(LISTV(1)))
        IF (abs(AREA) < 1.E-12) cycle
*
        NOF = NOF + 1
        LSTUF(NOF) = II
      end do
*
      write(80,"(a5,2i8)") 'CELLS', NOF, 4*NOF
      do II=1,NOF
        IFACE = LSTUF(II)
*
        call VEONF(IFACE,
     +  KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE,
     +  NOFVE,LISTV,LISTE)
*
        write(80,"(i1)",advance='no') NOFVE
*
        do JJ=1,NOFVE
          write(80,"(i8)",advance='no') LISTV(JJ)-1
        end do
        write(80,*)
      end do
*
*  ---  write cell_Types 
*
      write(80,"(a10,i8)")  'CELL_TYPES', NOF
      do II=1,NOF
        write(80,"(i1)") 5
      enddo
*
      close(80)
      return
      end