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

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

杉原厚吉先生の公開している2次元Delaunay法コードを用いて、2次元の点群を与えてDelaunay三角形分割を作成します。

2次元Delaunay

プログラムsiv2.fをコンパイルし、点群入力データdata.datを作成し実行すると上のようにうまくDelaunay分割が行われました。

しかし、一直線上に点が3点以上並ぶ場合は失敗しました。下図のように、三角形に接続していない点ができてしまいます。

2次元Delaunay失敗

これは、内部の計算が破綻した際に記号摂動法?で回避した結果と思われます。杉原先生のコードは計算が破綻しても決して落ちないロバストな設計となっています。

以下に、三角形分割の結果をVTK形式に出力するサブルーチンを示します。これをsiv2.fに追加するとParaView等で簡単に結果を見ることができます。

↑

DRDELVTK†

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
添付ファイル: filesugihara_del2d_２.png 3764件 [詳細] filesugihara_del2d_1.png 3561件 [詳細]
Last-modified: 2013-08-16 (金) 09:28:23