2004/09/10 山下 宏
１．このプログラムについて
このプログラムは二次元平面上において任意の形状をもつ物体に対して 有限要素法（Finite Element Method・・・FEM）で解析するための前処理を 行うものです。 対象となる物体の輪郭を入力してやれば、物体を自動的に三角形の要素 に分割し、要素番号、節点番号のデータを出力します。 特徴としては節点の密度分布に対応でき、ゆがみの少ないなめらかな三 角形データを作ることが出来ます。 16万個程度の節点を生成することが出来ます。 実際に解析するプログラムは別途、用意する必要があります。 １．はじめに 本プログラムはフリーソフトウェアとして公開致します。 ただし、本プログラムを使用したことによって生じた損害を補償する義務 は作者はまったく負いません。 ２．圧縮ファイルの内容 以下の６つのファイルが初期圧縮ファイルの内容です。 VORO95 .EXE VORO .HLP SAMPLE .DAT SUPANA .DAT VORO .DOC VORO_SRC.LZH VORO_SRC.LZH にソースリストが圧縮されて入っています。 かなり時間がかかるソフトなのでCPUは最低でも486以上のものが精神衛生上 必要となるでしょう。 プログラムは Visual C++5.0 でコンパイルされています。 製作機種はDOS/V (K6-III 450MHz)を使用しました。 ３．コピー、転載、移植、改良について コピー、転載は自由に行って構いません。 ただし、圧縮ファイルの中身を変更せずに配布して下さい。 特にSAMPLE.DAT は必ず同時に配布して下さい。 雑誌のCD-ROM等へ載せる場合はご一報頂けるとうれしいです。 このプログラムはソースリストも同時に配布しています。 移植、改良は自由です。また、移植、改良したものを再配布する場合は、 配付先を私までご一報下さい。 なお、最新版は http://www.yss-aya.com/voro.html に登録する予定です。 ４．おわりに ご使用にあたって、感想、バグ、要望等がございましたら yss あっと bd.mbn.or.jp （山下 宏）("あっと"を "@"に変換) まで御気軽にメールを頂ければ幸いです。２．使用方法
１．サンプルデータの使い方 それでは、まず SAMPLE.DAT を使って使い方の説明をします。 トップメニューの「ファイル」から「座標データの読み込み」を選んで下さい。 ファイル名を選択するダイアログボックスが開きますので "SAMPLE.DAT" を選ん で下さい。すると、図形の輪郭データを読み込み、SAMPLE.VOR というファイル を出力してもいいかと聞いてきますので「ＯＫ」を押して下さい。 コンピュータの速度にもよりますが Pentium の 90MHz 程度ならば１秒程で作成 が終わりデータが出力されるはずです。 細かく要素が生成されている所は拡大して見ることが出来ます。 拡大させたい場所を左クリックして下さい。最大65536倍まで拡大出来ます。 縮小は右クリックです。 もう少し、全体的に細かく切りたいな、と思った場合は座標データファイルを 修正します。 ここではトップメニューの「ファイル」「座標データの修正」を選んで下さい。 "メモ帳"が起動され、"SAMPLE.DAT"の内容が表示されます。 そこで基本密度半径の指定をしている場所を見つけ、 basic_r = 15.000 と、なっているのを basic_r = 8.000 に書き直します。直したら メモ帳のトップメニューの「ファイル(F)」「上書き保存(S)」で変更内容をセーブします。 再びこのプログラムのトップメニューの「ファイル」「座標データの読み込み」 から "SAMPLE.DAT" を選んで下さい。 今度は４秒程で（Pentium.90MHz)より全体的に要素が細かく切られたデータが得られます。 以上で、基本的な操作は終了です。 "SUPANA.DAT" はスパナのデータです。同様に選択すれば、以下、右の穴を細かく分割 したスパナが出力されます。 ２．その他のコマンドについて ヴォロノイ図の表示 ・・・ 水色でヴォロノイ図を上書きします。再描画されません。 三角形のヒストグラム・・・ 作られた三角形の要素の歪み具合を示すものです。 三角形の最小角が60度に近いほど（正三角形に近いほど）歪みが少ないと考え、 最小角ごとの要素の数を出力ファイルに書き出します。 ３度づつだと54〜58、57〜59度、というように３度単位でまとめて出力します。 通常は「なし」でご使用下さい。 ３．実際の入力データ作成 実際に自分でデータを作成する場合には SAMPLE.DAT を変更して作成して下さい。 ファイル名の拡張子は".DAT"でなくとも構いません。が、".DAT"にした方が便利です。 最低限、設定すべきデータは以下の通りです。 1. 物体の輪郭データ 2. 基本となる半径の大きさ 3. 要素を細かく切る点の指定 4. 点生成時のサーチ間隔の逆数（2.0〜5.0：整数にして下さい） 5. 分割したい点の場所 全てのデータは浮動小数であることを自明とするために 455.0 のように「.0」付きで 指定してやった方が安心です。 1. の輪郭データはｘｙ平面上に座標データの形で入力します。 その際、データの最後は必ず始点と同じ座標を与えて、そこで図形が「閉じる」ことを 明確にしてやる必要があります。 図形内部に空洞がある場合は、さらに続けてデータを入力して下さい。 ただし、図形同士が交差する、といった通常はありえない設定に対してのエラー チェックはしていません。 また、全ての図形は連絡しているものとします。図形の中に図形を定義した場合は その内部は空洞と見なされます。その空洞の中にさらに図形を定義することは意味が ないと思われますが、その場合、その内部は実際にあるものとして認識されます。 なお閉じた図形の最大数は256となっています。 また、ｘｙ座標データは1024個を超えて入力することは出来ません。 2. の基本となる半径の大きさは、全体の大きさを100×100とした場合、 10を指定すればノードは約100生成され要素は200生成されます。 あまり小さくしすぎると、処理が遅くなります。 基本半径には最初は大きな値を与え、徐々に小さくしていくのがいいでしょう。 なお、大きすぎると、内部の輪郭に対して、正しく要素を生成しない場合が あります。 実際にデータを使用する前にディスプレイ上で正しく要素が生成されてる事を 確認してから使用することを勧めます。 3. は要素を細かく切りたい場所を指定するものです。別に指定しなくても 構いません。与え方は１つの指定点につき１行で指定します。 順番は ｘ ｙ 最小密度（半径） それが及ぼす半径 の順で指定します。 最小密度というのは（ｘ，ｙ）地点における要素の半径（外接円）を表しています。 これが小さいほど、そこの場所の要素は細かくなります。 それにつづく半径も単位としては同じです。 これを小さくすれば（ｘ，ｙ）にごく近い付近にしか影響を与えません。 4. の、点生成時の間隔は通常はあまり変更する必要はありません。 本プログラムでは点を生成する場合に画面の下から上に向かってある間隔で座標を 調べ、そこに点をうっても大丈夫かを判定しています。そのため、間隔が狭いほど正 確に点を打っていきます。この値は3.0ならばその位置の半径の3分の1の間隔でサーチ することを意味します。間隔を 2.0 にすると密度が高い場所では、かなり荒っぽく点 が打たれているのが分かります。 通常は 2.0〜5.0 で大丈夫だと思いますが、ノードの数を多くする場合などには 5.0 だと時間がかかるので 2.0 にするといいかもしれません。 5. の分割したい場所はｘｙ座標データの形で入力します。指定しなくても構いません。 この点を指定すると、この場所には必ずノードが生成されます。 ただし、必ず図形の内部（境界線上を含む）にして下さい。それ以外の場所では動作は保証できません。３．出力データの説明
出力されるファイルは、入力ファイルの拡張子を ".VOR" に変更した物です。 拡張子がない場合には ".VOR" を付けます。この出力ファイル名は変更出来ません。 出力ファイルは入力ファイルと同じディレクトリに置かれます。 出力されるデータは以下の順番です。 基本的にそれぞれのデータは空白で区切られています。 1. ノードのｘｙ座標 2. 要素（メッシュ）の接続する３点 3. 輪郭指定点のノード番号と座標 4. 輪郭上の点と座標 5. 分割指定点のノード番号と座標 1. では最初に全ノードの個数を出し、続いて「ノード番号 ｘ座標 ｙ座標」の 順番で出力します。 2. では最初に全要素の個数を出し、続いて３点のノードの番号をだします。 要素データは反時計回りになるようにノードは並び替えられています。 3. では輪郭を指定したノードの番号とｘｙ座標を出力します。 閉じた図形ごとに１行開けて出力します。 4. でも輪郭上の点を出力し、閉じた図形ごとに１行開けて出力します。 5. では分割指定点のノード番号と座標を出力します。４．プログラムの説明
このプログラムの最大の「売り」は点の打ち方にあります。 概念としては平面において「任意の半径の円を敷き詰める」という考えで行って います。すなわち、節点の間隔が狭いところでは小さな円を敷き詰め、間隔の粗い 場所では大きな半径の円を敷き詰めます。これによって、点の細かいところから粗 いところまで比較的なめらかに点を打つことが可能となっています。 まず、全体の密度（半径）分布が与えられ、それを満たすように点を配置してい きます。 最初に輪郭指定点には必ず点を打ちます。そして、次に輪郭線上の点に対して、 その位置における密度にあった間隔で節点を打ちます。この場合、輪郭が鋭角にな っていて節点の密度が関数値を超える場合でも生成します。 輪郭指定点から次の輪郭指定点まで点を打っていき、次の点にぶつかると、中点 をとって、あまり歪まないように補正しています。 そして、形状内部への節点の配置方法は、形状の一番下からある間隔でサーチし ていき、ある位置に点を打てるかどうかは、その位置からもっとも近い既存の点ま での距離と、その位置での半径を比較して距離が半径より大きいならば点を打つ、 というルールです。 実際のプログラムでは点を打つと同時にヴォロノイ図をも形成していくため、 最近点を求める計算は苦になりません。 本来ならばバケット法を用いて計算途中で消えていくヴォロノイ点や辺の数を減 らすべきなのですが、点を打つ方法との兼ね合いのため、そうしていません。 サーチする点の間隔を細かくすると、より正確に点を打つのですが、速度が問題 となるため実際は、その場所の密度の２分の１から５分の１の間隔でサーチするよ うにしています。（このパラメータがスパンです） ここで、ヴォロノイ図という言葉が出てきましたが、この図はある点が幾つか与 えられた場合、その点から最も近い領域を示す図で、その境界は点同士の垂直二等 分線で区切られた多角形となります。 なお、三角形要素形成はヴォロノイ図を形成しdelaunay法（ヴォロノイ図におい て境界線をもつ点同士を結ぶ）を用いてゆがみの少ないように節点を線で結んでい きます。 この際、もっともやっかいなのが数値計算誤差による図形判定処理アルゴリズム の破綻と４個以上の母点が同一円周上にある場合、いわゆる退化が生じた場合の処 理です。 これらの現象が生じるとまず、プログラムが無限ループに陥るなどの症状を示し、 また終了した場合でも、どのヴォロノイ点にも属さないヴォロノイ線が残り、誤っ た節点を結んでしまい、正常な三角形要素が作れないなどの症状が現れ、システム 全体が非常に不安定な状態になります。 初期段階においては、まず退化を起こさないような節点を作成し、なおかつ、処 理終了後、正しく要素が形成出来ているかを、ディスプレイ上で確認して、ごまか してデータを出力していました。 しかし、実用にならないので参考文献の杉原厚吉氏のアルゴリズムを用いました。 これはヴォロノイ図を位相構造のみを正確に計算することによって作るアルゴリズ ムで、このアルゴリズムにより数値計算誤差の影響を減らし、誤動作が極めて少な くなるようになりました。 なお、本プログラムにおいては杉原氏のアルゴリズムは完全には再現されていま せんので、誤動作する危険があります。 （現時点においてはまだ誤動作は見られません） ヴォロノイ図等の詳しい説明は参考文献をご覧下さい。５．制限事項
・生成可能なノード（節点）の限界数は16万個程です。 ただし、それなりのメモリと計算時間がかかります。 例として supana.dat において基本密度半径を1.0に、密度データの最小密度を0.1にした場合には ノード数 = 36017 計算時間 = 412秒 使用メモリ = 13MB （メモリは最低でも32MB以上必要でしょう。スワップしまくる） 出力ファイルサイズ = 2MB 以上、Pentium90MHz（RAM72MB）のマシンにおいて計測しました。 ノード数が少ない時は、それなりのメモリしか使用しませんので問題ありません。 （Win95のメモリ管理が優秀なため） ・大きな入力ファイルを読み込もうとすると一般保護違反で落ちる場合があります。通常はVORO用のデータファイルではない物を読み込もうとした場合です。 ・図形を拡大しすぎると、一部の線が消える場合があります。（ int の範囲を超えるため） ・計算中は入力を全く受け付けません。 ・Voronoi図を表示しようとした場合にハングアップする場合があります。 Line文で負の大きな値が入った場合に起こります。原因不明。クリッピングミスか？ ディスプレイドライバによるのかもしれません。６．コマンドラインでの指定
・コマンドラインで voro95 sample.dat のようにファイル名を指定して実行するとエラーがない場合は 画面を表示することなくデータを作成して終了します。７．参考文献
１．矢川 元基・吉村 忍 有限要素法へのニューロ・ファジィの応用 シミュレーション 第13巻第1号 (1994) ２．伊理 正夫監修／腰塚 武志編集 計算幾何学と地質情報処理 第２版 共立出版株式会社 (1993) ３．杉原 厚吉 幾何アルゴリズムの数値的破綻とその対策 応用数理 Vol.1 No.4 DEC. (1991)８．履歴
Version 1.00 1995年3月8日 暫定公開バージョン。 取りあえずまともに動くようなのでフリーソフトとして公開する事にした。 と、いうよりも卒論の締切が迫ってきたので無理矢理形にした(^_^;) Version 1.01 1995年3月17日 386専用コードで生成されていたのを8086コードに変換した。 Version 2.00 1996年7月10日 9801用のプログラムをWindows3.1用に書き換える。 かなり使いやすくなった、と自分では思うが・・・。しかし、Windowsの ソフトを作るのはとっても面倒であーる。MFCを使えば楽になるのかな？ Version 2.03 1996年7月12日 三角形のひずみ具合を示すヒストグラムオプションを付ける。 拡大率を65536までに。 密度指定を逆にしたときに出る０割エラーを修正。 図形が横長の時の無駄なサーチの打ち切り。 Version 3.00 1996年7月12日 Windows95専用版。32ビット化にともない、可能なノード（節点）の数が １万を超える。また計算速度も95上で16ビット版の1.9倍も早く動作した。 （huge array を使っていたせいかな？） Version 3.01 1996年7月20日 全ての float を double にしたことにより、より細かいメッシュを切った 時の、要素作成ミスを減らす。 Version 3.02 1998年2月27日 ｙ軸に垂直な境界線上で、ノードがまれに時計回りで生成されてしまっていたのを訂正。 河原由夫さん、ご指摘ありがとうございました。 Version 3.10 1998年4月1日 分割指定点を追加。沙培正さん、要望ありがとうございました。 Version 3.11 1998年6月9日 三角形の歪み具合を算出する部分で保護違反で落ちることがあるのを訂正。 Version 3.12 1998年10月26日 出力データが全て時計回りだったのを訂正(河原由夫さん指摘)。なぜ作った時に気がつかなかったのだろう・・・。 Version 3.13 2000年01月14日 *.dat データファイル中に空白だけの行があっても読みこめるように。 Version 3.14 2000年01月15日 たまに時計回りの三角形ができるのを修正？したつもり。 Version 3.15 2000年12月10日 元座標データが重複していた場合に警告するように。 Version 3.16 2004年09月10日 擬似的にコマンドラインでファイル名を指定した場合に自動でデータを作成して終了するように。９．謝辞
この研究を行うにあたって御指導頂いた東北大学地球工学科 千田 佶教授ならび に洪 承燮助教授に感謝致します。 また、研究のみならず、多々の相談にものっていただいた新堀雄一助教授に感謝 致します。 また、資源システム工学講座(７講座）の皆さんには、研究室での生活を楽しく、 有意義なものとさせて頂いた事を紙面をお借りして感謝いたします。 皆様方の今後の御活躍と御健康を祈りつつ、謝辞と代えさせていただきます。 −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− なお、このヘルプファイルの作成には、 「さるヘル」・・・本名「サルでも書けるヘルプファイル」 梅木泰宏氏作フリーソフト 登録場所:FGALAM #13, JBC03637 を使わせていただきました。とても使いやすいソフトなのでお勧めです。