*[[2次元Delaunay法(杉原厚吉先生)]] [#a59dbce3]

[[杉原厚吉:http://ja.wikipedia.org/wiki/%E6%9D%89%E5%8E%9F%E5%8E%9A%E5%90%89]]先生の公開している[[2次元Delaunay法コード:http://home.mims.meiji.ac.jp/~sugihara/kakekomi/iwanamisoft/iwanamisoft.html]]を用いて、2次元の点群を与えてDelaunay三角形分割を作成します。

#ref(sugihara_del2d_1.png,center,nowrap,80%,2次元Delaunay);

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

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

#ref(sugihara_del2d_2.png,center,nowrap,70%,2次元Delaunay失敗);

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

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

**DRDELVTK [#y88ae58d]
#code(fortran){{
      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

}}

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS