JPH08315183A

JPH08315183A - 自動メッシュ生成方法及びシステム

Info

Publication number: JPH08315183A
Application number: JP7117244A
Authority: JP
Inventors: Kenji Shimada; 憲司嶋田
Original assignee: IBM Japan Ltd
Current assignee: IBM Japan Ltd
Priority date: 1995-05-16
Filing date: 1995-05-16
Publication date: 1996-11-29
Anticipated expiration: 2014-04-12
Also published as: JP2881389B2

Abstract

(57)【要約】【目的】非多様体データ・モデル上で、物理モデルに
よる自動メッシュ分割を適用するためのシステム及び方
法を提供すること。【構成】非多様体データ・モデルを用意し、所定領域
内に複数の楕円バブルを初期的に配置する。この楕円バ
ブルは、与えられたベクトル場のそれぞれのベクトルに
対応した大きさ及び方向性を有している。このベクトル
場は、流体の流れなどといったある自然現象を示すもの
である。楕円バブルを定義されたバブル間力により移動
させるとともに、バブルの個数の適応的に制御すること
により、所定領域内に過不足なく配置する。このように
楕円バブルが安定した時点で、Delaunay法を用いて、バ
ブルの中心点を結ぶことによって、自然現象を示すベク
トル場に対応した三角形に分割することができる。

Description

【発明の詳細な説明】【０００１】【産業上の利用分野】この発明の目的は、有限要素法、
ＣＡＤシステム、コンピュータ・グラフィックスなどの
技術分野において適用される自動三角要素分割技法に関
するものである。【０００２】【従来の技術】近年、形状の美しさと機能の向上を目的
として、自動車や家庭電化製品を始めとする製品形状に
複雑な自由曲面が多く使われるようになってきた。一般
に、これらの曲面は設計時には Bezier、Ｂ−スプライ
ン、NURBS などのパラメトリック曲面としてＣＡＤシス
テム内で表現されるが、他方、解析のために有限要素法
などを用いて製品性能を確認する際には別の形状表現、
すなわちメッシュが必要となる。したがって設計と解析
の工程を有機的に結ぶためには、ＣＡＤにおける曲面形
状から三角形メッシュへの効率良い自動変換の方法が求
められている。【０００３】従来より、この技術分野の刊行物として次
のようなものが知られている。【０００４】すなわち、特開昭３−１６３６６９号公報
は、ＣＡＤの分野で図形解析などに用いられるメッシュ
発生方法に関し、特に、２次元又は３次元領域に点を逐
次追加して三角形または四面体要素でデローネイ（Dela
unay）分割を行うメッシュ発生方法に関し、２次元又は
３次元領域の領域を予め規則的パケット要素に分割し、
新たなメッシュ点を追加した際、接続が組み換えられる
要素群の検出を、その点が所属するパケットに外接円ま
たは外接球が登録されている要素に限って行うことによ
り、組み換え要素検出を高速に行うようにすることを開
示する。【０００５】特開平１−２８６０８１号公報は、多次元
オブジェクトの多角形ドメインの関する情報を入力し
て、この情報を処理し、上記処理の結果に従って表示を
発生するタイプの、予測された物理現象の多重多角形表
示を発生する方法において、直線の辺もしくは曲線部分
である境界を有し、各々が少なくとも１点で上記ドメイ
ンの境界と接する、複数の第１の領域に、上記ドメイン
を分割し、任意の曲線部分を、該曲線部分の両終点を接
続する少なくとも一本の直線で置換して、上記第１の領
域を調整し、これによって該調整された第１の領域のす
べての辺が線分で囲まれるようにし、共通の境界によっ
て囲まれた、上記調整された第１の領域の対を有限数の
クラスに分類し、該各々のクラスに適した規則に従っ
て、上記第１の領域を第２の領域に分割して、該第２の
領域を上記発生の表示に使用することを開示する。【０００６】上記刊行物に示されているものも含め、こ
れまで、三角要素分割問題は計算力学の分野において最
も多く研究・開発されてきた。実際、２次元形状につい
ては多くの方法が提案されているが、３次元形状につい
ては必ずしも十分に研究されていないのが現状である。【０００７】また、従来の方法は主として２次元・３次
元形状領域に特化しており、非多様体形状に対して一般
的に適応できるものではない。これはいくつかの市販の
ＣＡＤシステムにおいて３次元非多様体形状が表現でき
るようになったのが比較的最近のことだからであろうと
思われる。【０００８】さて、従来の典型的な三角要素分割法を、
K. Ho-Le, "Finite element mesh generation methods:
a review and classification," Computer Aided Des
ign,Vol.20, No.1, 1988に基づき、次の４グループに分
けて説明する。【０００９】部分領域分割法与えられた領域を凸領域や穴のない領域などの単純で扱
いやすい部分領域に切り分けて、次にそれぞれの部分領
域を三角要素分割する。部分領域の分割は、規則正しい
格子を写像関数でマッピングしたり、簡単なルールによ
って内部ノードを加えたりして行う。従来実現されてい
る方法では、部分領域への分割を人手に頼っていること
が多く、完全には自動化されていないものがほとんどで
あった。また、それぞれの部分領域に指定されたノード
間距離分布に忠実にノードを配置する方法がなかった。【００１０】階層的空間分割法２次元領域を４分木で、３次元領域を８分木で階層的に
再分割してゆく方法．滑らかな領域境界を表すために、
境界近傍の要素として一部が欠けた正方形や立方体を扱
えるように拡張する方法も提案されている。最大の問題
点は分割の深さによって要素の大きさを離散的にしか制
御できないことであろう。また、分割軸に対する形状の
置きかたが分割結果に大きく影響することや、解析など
で重要となる形状の角や表面に歪みの大きい要素が数多
く生成されるなどの問題点も指摘されている。【００１１】再帰的二分割法三角要素として使える形状が得られるまで、領域を再帰
的に二分割してゆく方法．部分領域分割法と異なり、要
素のレベルまで二分割だけで領域を切り分けていく。こ
の際に必要となる形状モデリングのデータ操作などの議
論は多く行われているが、ノード間距離の連続的な制御
や要素の歪みを少なくすることなどの基本的要求に対し
てはまだ十分に考慮されていない。【００１２】ノード結合法まず領域の境界と内部にノードを配置し、次にこれらを
結んでメッシュを構成する方法。いったんノードの集合
が与えられると、Delaunayの方法を用いてこれを結合し
て歪みの少ないメッシュを生成できる。この方法は計算
幾何学の分野で長く研究されてきたもので、２次元・３
次元空間に散らばった点を接続して正三角形や正四面体
に近いメッシュを生成するための効率的かつ安定なアル
ゴリズムである。（例えば、伊理，計算幾何学と地理情
報処理，bit 別冊，共立出版，1989を参照）。残された
問題はいかにして指定されたノード間距離分布を満たし
歪みの少ない要素を生成するようなノード配置を得るか
である。特に３次元の場合は難しい。【００１３】ところで、有限要素法などを用いた解析に
おいて、計算時間や記憶空間を増やさずに解精度を改善
するには、領域全体にわたりノード間距離（すなわち三
角形の辺の長さ）を連続的な分布関数にしたがって変化
させなければならない。この際、解析解の値が急に変化
する領域では細かな三角要素（三角要素という用語は、
平面における三角形要素と、空間における四面体要素の
総称である）を、それ以外では粗い要素を使うことが大
切となる。【００１４】また、与えられたトリミング曲線上や、特
に指定された点や曲線上にノードが正確に置かれている
ことが必要である。また、解析においてこれらの曲線上
で境界条件を与えることができるように、各曲線の近似
形状が三角形の辺の列としてメッシュ内に表現されてい
なくてはならない。【００１５】【発明が解決しようとする課題】一般に、メッシュを構
成する要素は、正三角形や正四面体等の等方的なメッシ
ュを用いることが好ましいとされてきた。これは、著し
く歪んだメッシュは解析の精度を低下させるからであ
る。しかしながら、いくつかの解析問題においては、等
方性メッシュではなく、ある自然現象と関連性を有する
ベクトル場の方向に対応した細長く伸びた非等方性メッ
シュを用いる方が、解の精度や収束性の改善という点
で、好ましい場合がある。例として、自動車の衝突解析
や板成形のシミュレーションなどの大変形の弾塑性変形
解析があげられる。衝突解析においては、梁部材に力が
流れる方向に細長く伸びたメッシュを用いるのがよい。
すなわち、力の流れる方向というベクトル場に沿って、
非等方的なメッシュを用いることが好ましい。また、流
体解析でも流線の方向というベクトル場に沿って細長く
伸びたメッシュを用いる方が解析の効率がよい。【００１６】図１８は、流線のベクトル場の一例を示す
図である。この例におけるベクトル場は、自動車の風洞
実験における風の流れであり、この風の流れに応じてメ
ッシュを非等方的に形成することが好ましい。また、構
造解析を行う場合、ベクトル場は力の流れであり、曲面
の近似においては、最大曲率ベクトルである。【００１７】このように、ある自然現象に関係するベク
トル場に応じて、メッシュを非等方性にする理由は、自
然現象を表すベクトルと解となる物理量の変化に一定の
関係があることに基づく。すなわち、「力の流れの方
向」、「板が曲げられる方向」或いは、「流線の方向」
などといった自然現象を表すベクトルの方向と同一方向
に対しては、解となる物理量の変化が少なく、これと直
交する方向は、物理量の変化が大きいという先験的知識
に基づいている。このような関係ゆえに、自然現象を示
すベクトルの方向と同一方向については、物理量の変化
が少ないためメッシュを大きくしても解の精度を低下さ
せることがないので構わないが、これと直交する方向に
ついては、物理量の変化が大きいためメッシュを小さく
することが好ましいわけである。【００１８】そこで、この発明の目的は、自然現象を示
すベクトル場を考慮して、解の精度が高く処理速度の速
い自動メッシュ分割を適用するためのシステム及び方法
を提供することにある。【００１９】この発明の他の目的は、非多様体データ・
モデルに複数の楕円バブルを配置し、該複数の楕円バブ
ルの間での定義された引力・斥力に基づき力学的な方程
式を解き、楕円バブルの中心点を結ぶことにより、自動
的に最適な三角要素分割を実行するシステム及び方法を
提供することにある。【００２０】【課題を解決するための手段】上記目的は、次に示すよ
うな本発明のステップ、及びそれらのステップを実現す
るシステムによって達成される。【００２１】本発明は、予め用意された非多様体データ
・モデルに対して、以下に示すステップ１乃至ステップ
５の手順により最適な数のノードを最適な位置に発生す
る。また、以下のステップ６に示したようにこのノード
発生方法をワイヤーフレーム（稜線）、サーフェス（曲
面）、ソリッド（空間領域）の順に適用し、さらに三角
要素に結んでいくことにより、３次元非多様体の自動三
角要素分割を実現するものである。なお、図１は、本発
明の概要を示す図である。【００２２】ステップ１．自然現象を示すベクトル場の
指定メッシュの特性を表す基本となるベクトル場を指定す
る。すなわち、要素の大きさを指定するスカラー場ｄ及
びメッシュの方向性とアスペクト比を指定するベクトル
場ｖを指定する（図１（ａ））。【００２３】ステップ２．楕円バブルの定義所望の三角要素が得られるような理想的なメッシュ・ノ
ード配置を求めるために、各ノードを上記のスカラー場
ｄとベクトル場ｖに応じて大きさと方向が決められる楕
円（以下、楕円バブルという）として表現する（図１
（ｂ））。【００２４】ステップ３．楕円バブル間力を定義楕円バブルを形状領域内に最密に充填するために、隣接
する二つの楕円バブル間に働く力を定義する。すなわ
ち、楕円バブル間には分子間力に似た力の場をバブル間
距離の関数として定義する。二つの楕円バブルは、接し
ているとき等、二つの楕円バブルの距離がある所定の安
定距離内にある場合が安定状態であり、楕円バブル間力
は働かない。そして、この距離よりも近づくと斥力が働
き、遠ざかると引力が働くようにする。好ましい楕円バ
ブル間力は、例えば三次曲線の一部分などにより定義さ
れる。【００２５】ステップ４．動力学シミュレーションによ
る楕円バブル間力の均衡位置の算出本ステップの動力学シミュレーションと以下のステップ
の適応的楕円バブルの個数制御により、適切な楕円バブ
ルの個数と全ての楕円バブル間力が均衡して安定する配
置を求める。この状態では、楕円バブルが最密充填され
た状態になる。このような状態にするために、上記の楕
円バブル間力に加えて、楕円バブルの中心点に質量と楕
円バブルの速度に比例する粘性抵抗を与えて古典的ニュ
ートン力学の運動方程式をたてる。ある初期位置から始
めてこの運動方程式をRunge-Kutta法などの数値積分に
よって解くと、時間の経過とともに楕円バブルの配置が
力の均衡状態に近づいてゆくので、十分に近づいたとこ
ろで数値積分をうちきり最終的な解とする。ただし、ワ
イヤー・フレーム上とサーフェス上の楕円バブルはそれ
らの幾何要素上から楕円バブルが外に出ないように移動
させなければならない。【００２６】ステップ５．楕円バブル個数の適応制御上記の動力学シミュレーションの最中に楕円バブル同士
の重なり状態を計算して、この重なりの程度を重なり度
として定義する。そして、この重なり度に応じて、楕円
バブルの密度が大きすぎる場所では適宜、楕円バブルを
破壊して個数を減らし、小さすぎる場所では分裂させて
個数を増やす。このメカニズムによってある領域に対し
て適切な数の楕円バブルを過不足なく充填することがで
きる（図１（ｃ））。【００２７】ステップ６．非多様体形状上への楕円バブ
ルの配置と三角要素分割最密充填された楕円バブルの中心点をメッシュ・ノード
とし、これらのメッシュ方向性とアスペクト比とを考慮
した「デローニ三角形の判定条件」を用いて、各メッシ
ュ・ノードを自然現象を示すベクトル場の方向に細長く
なるように結んで三角メッシュを形成する（図１
（ｄ））。幾何要素に対して上述の楕円バブルの発生方
法が用意されているとき、３次元非多様体形状の三角要
素分割は以下の手順で行うことができる。【００２８】１．全ての頂点に楕円バブルを配置する。２．全ての稜線上に楕円バブルを最密充填する。３．稜線上の楕円バブルの中心点を順に結んで線分要素
に分割する。４．全ての面上の閉じた領域に楕円バブルを最密充填す
る。５．面上の閉じた領域内の楕円バブルの中心点を、後に
詳述する変更されたDelaunayの方法などで結んで三角形
要素に分割する。６．全ての空間内の閉じた領域に楕円バブルを最密充填
する。７．空間内の閉じた領域内の楕円バブルの中心点を変更
されたDelaunayの方法などで結んで四面体要素に分割す
る。【００２９】尚、特に非多様体モデルにあっては、上記
楕円バブルの充填順序を上述の順序で行うことが望まし
い。さらに、本発明にあっては、非多様体モデルの通常
の稜線と、次元の異なるオブジェクト間の境界をなす稜
線を区別せず同様に楕円バブル充填を行うことが望まし
い。また、最初に配置される初期楕円バブルは、規則的
格子の位置に置いてもよいし、階層的な領域分割（２分
木、４分木、８分木など）を用いてもよい。【００３０】【実施例】以下図面に基づき、本発明の好ましい実施例
について説明する。【００３１】図２は、本発明の自動三角要素分割システ
ムを実現するための処理モジュールのブロック図であ
る。この実施例では、図２に示す各々の処理モジュール
のブロックは、ＡＩＸ（ＩＢＭの商標）の下で動作する
ワークステーション上で、Ｃ言語で作成された個別のプ
ログラム、ルーチンまたは関数である。しかし、これら
個々のモジュールを、ディスクリートなハードウェアで
実現してもよいし、複数台のワークステーションをＬＡ
Ｎなどで相互接続した分散処理環境で、個別のワークス
テーションに別の処理モジュールをロードし、分散並列
処理を行うようにしてもよい。【００３２】図２において、形状入力ブロック１は、メ
ッシュ分割すべき入力形状を定義し、要素分割に際して
必要な情報を提供するもので、メッシュ分割すべき入力
形状を定義する非多様体入力形状定義手段１０２と、こ
の形状から、これらを構成する頂点、稜線、面、空間な
どの位相要素を得る入力形状位相情報抽出手段１０４
と、これらの位相要素の位置と形状などの幾何情報を得
る入力形状幾何情報抽出手段１０６とから構成される。【００３３】形状入力ブロック１は、スカラー場定義手
段１０７及びベクトル場抽出手段１０８を有している。
スカラー場定義手段１０７により、要素の大きさを指定
するスカラー場ｄが定義される。またベクトル場定義手
段１０８により、メッシュの方向性とアスペクト比とを
指定するベクトル場ｖが定義される。このスカラー場ｄ
及びベクトル場ｖはそれぞれ、２次元ではｄ（ｘ、
ｙ）、ｖ（ｘ、ｙ）で表され、３次元ではｄ（ｘ、ｙ、
ｚ）、ｖ（ｘ、ｙ、ｚ）で表される。なお、このベクト
ル場は、上述のように力の流れの方向や流線の方向等の
ある自然現象の方向を表すものである。【００３４】要素分割ブロック２は、形状入力ブロック
１から順次、必要な位相的・幾何的な情報を得て、メッ
シュ・ノードを配置するとともにメッシュ要素を発生
し、さらに生成されたメッシュ出力ブロック３に順次出
力する。【００３５】より具体的には、要素分割ブロック２は、
形状入力ブロック１からの情報に基づき頂点のノードに
楕円バブル（詳細は後述する）を配置する頂点ノード配
置手段２０２と、形状入力ブロック１からの情報に基づ
き稜線上に楕円バブルを配置する稜線ノード配置手段２
０４と、配置された頂点及び稜線から線要素を発生する
線要素発生手段２０６と、形状入力ブロック１からの情
報に基づき面ノード上に楕円バブルを配置する面ノード
配置手段２０８と、配置された面ノードから２次元のDe
launayの方法に基づき三角形要素を発生する三角形要素
発生手段２１０と、形状入力ブロック１からの情報に基
づき空間ノード上に楕円バブルを配置する空間ノード発
生手段２１２と、配置された空間ノードから３次元のDe
launayの方法に基づき四面体要素を発生する四面体要素
発生手段２１４からなる。特に非多様体データ構造の三
角形分割を適切に行うために、要素分割ブロック２にお
ける各手段２０２〜２１４は、図１に示す上から下に順
次実行する。【００３６】メッシュ出力ブロック３は、要素分割ブロ
ック２から出力されるノードと、メッシュ要素の位相的
・幾何的情報を記憶しておき、必要に応じて要素分割ブ
ロック２の処理過程に、保持していた値をフィードバッ
クする。このため、メッシュ出力ブロック３は、ノード
とメッシュ要素の位相的・幾何的情報を記憶しておくメ
ッシュ定義手段３１０と、メッシュ定義手段３１０に位
相的情報を書き込むメッシュ位相情報生成手段３０２
と、メッシュ定義手段３１０に幾何的情報を書き込むメ
ッシュ位相情報生成手段３０４と、メッシュ定義手段３
１０から情報を得て要素分割ブロック２に位相的情報を
フィードバックするメッシュ位相情報抽出手段３０６
と、メッシュ定義手段３１０から情報を得て要素分割ブ
ロック２に幾何的情報をフィードバックするメッシュ幾
何情報抽出手段３０８から構成される。【００３７】次に、図３及び図７を参照して、本発明の
処理ステップについて説明する。図３は、本発明に係る
処理の概要ステップを示すフローチャートである。図３
では、図２の形状入力ブロック１によって、図７（ａ）
に示すような非多様体入力形状が与えられているものと
する。尚、図７（ａ）において、３次元図形である四面
体の立体６００２に、稜線６００５を境界線として２次
元図形である三角形６００４が接し、さらに三角形６０
０４に、一次元図形である直線６００８が一点６００７
を共有して接している。このような異なる次元のオブジ
ェクトからなる複合図形は、非多様体データ構造を使用
しないと表現が困難であるかまたは実質的に不可能であ
る。また、立体６００２の１つの側面にはやはり一次元
の図形である「×」形状の図形が張り付いており、さら
には、立体６００２は、自身を分ける境界線６０１２を
内在している。このような形状を保持することも、非多
様体データ構造の１つの特徴である。【００３８】図３において、ステップ２００２では先
ず、図２の頂点ノード配置手段２０２と、稜線ノード配
置手段２０４によって、稜線上の頂点と、両端点及び内
部点と、稜線上に楕円バブルが配置される（図７の
（ｂ）及び（ｃ）を参照）。ここで、３次元の非多様体
データ構造において注意すべきなのは、点は稜線、面及
び空間の全てにおいて境界になり得るものであり、稜線
は、面及び空間の両方において境界になり得るものであ
り、面は、空間において境界になり得るものである、と
いうことである。以下の説明では、「境界」とは、この
意味で使用される。例えば、非多様体データ構造では、
空間中に、孤立した点が単独で存在し得るのであって、
この場合、その点自体が、周囲の空間との境界をなして
いる。【００３９】ステップ２００４では、後に詳述するよう
な力学的技法によって楕円バブルが移動され、必要に応
じて楕円バブルが破壊・分裂される。ステップ２００６
では、楕円バブル位置にメッシュ・ノードが配置され
る。ステップ２００７では、線要素発生手段２０６によ
って、楕円バブルの中心位置を繋ぐように線要素が発生
される。尚、ステップ２００２〜２００７は、図４のフ
ローチャートに関連して詳述される。【００４０】ステップ２００８では、図２の線要素発生
手段２０６によって、線要素を発生した後（図７
（ｄ））、面ノード配置手段２０８によって面の境界上
と面上に楕円バブルが配置され、ステップ２０１０では
上記ステップ２００４と同様の技法によって楕円バブル
が移動され、必要に応じて楕円バブルが破壊・分裂され
る。ステップ２０１２では、楕円バブル位置にメッシュ
・ノードが配置される（図７（ｅ））。すると、ステッ
プ２０１３では、三角形要素発生手段２１０により、後
に詳述する修正された２次元のDelaunayの方法を利用し
て、楕円バブルの中心位置が繋がれ、以って三角形要素
が発生される（図７（ｆ））。このとき注意すべきであ
るのは、非多様体データ構造においては、図７（ａ）の
「×」形状６０１０で示されるように、２次元の面上に
１次元の線が埋めこまれていたりすることである。従っ
て、修正された２次元のDelaunayの方法を適用する際に
は、このような埋めこまれた１次元の線に交差しないよ
うに、三角形要素を発生する必要があることに留意され
たい。ステップ２００８〜２０１３は、図５のフローチ
ャートに関連して詳述される。【００４１】ステップ２０１４では、図２の三角形要素
発生手段２１０によって三角形要素が発生された後（図
７（ｆ））、空間ノード配置手段２１２によって空間の
境界上と空間の内部に楕円バブルが配置され、ステップ
２０１６では上記ステップ２００４と同様の技法によっ
て楕円バブルが移動され、必要に応じて楕円バブルが破
壊・分裂される。ステップ２０１８では、楕円バブル位
置にメッシュ・ノードが配置され、ステップ２０１９で
は、修正された３次元のDelaunayの方法を利用して、そ
の空間内のメッシュ・ノードに関連して、四面体要素発
生手段２１４によって四面体要素が発生される（図７
（ｇ））。この場合にも、非多様体であるが故に、空間
内に埋めこまれた、空間よりも低次元の領域に交差しな
いように四面体要素を発生する必要があることに留意さ
れたい。ステップ２０１４〜２０１９は、図６のフロー
チャートに関連して詳述される。なお、変更されたDela
unayの方法については、後に詳述する。【００４２】このようにして結局、図７（ｈ）に示すメ
ッシュが生成される。【００４３】次に、図４以下を参照して、本発明に係る
メッシュ生成の詳細な処理について説明する。【００４４】図４を参照すると、稜線上の頂点及び稜線
上に楕円バブルを配置し、動力学的モデルに従い楕円バ
ブルを移動させ、必要に応じて楕円バブルを破壊・分裂
させるための処理のフローチャートが示されている。【００４５】まず、本発明で使用する楕円バブルについ
て説明する。座標（ｘ、ｙ、ｚ）に位置する楕円バブル
は、定義されたスカラー場ｄ≡ｄ（ｘ、ｙ、ｚ）に応じ
た大きさを有し、メッシュの方向性及びアスペクト比と
して指定されたベクトル場ｖ≡ｖ（ｘ、ｙ、ｚ）に応じ
た方向を有する楕円であって、質量ｍを有するような仮
想的な物体である。このスカラー場ｄ（ｘ、ｙ、ｚ）及
びベクトル場ｖ（ｘ、ｙ、ｚ）は、一般的には、定数で
はなく、空間の座標値に依存する関数として、それぞれ
スカラー場定義手段１０７、ベクトル場定義手段１０８
によって定義されている。例えば、建築物のＣＡＤデー
タのメッシュ生成を行う場合、応力が大きい箇所では、
比較的密なメッシュを生成させるように、ｄの値が小さ
くなる。ただし、スカラー場ｄというのは、楕円バブル
の近接度合い及び楕円バブル間力を示す指標として使用
されるものである。そして、ベクトル場ｖは、力の流れ
等のある自然現象を示している。従って、図８（ｃ）に
示すように、１つの楕円バブルは、他の楕円バブルと重
なり合わないようにも存在し得るし、互いに丁度接触す
るようにも存在し得るし、所定範囲以内にも近接し得
る。【００４６】楕円バブルは、基本となる円状のバブルの
大きさがスカラー場ｄによって決定され、さらにベクト
ル場ｖの方向にアスペクト比｜ｖ｜倍に引き伸ばすこと
によって、その大きさと形が定義される。【００４７】また、楕円バブル同志は、分子間力として
のファン・デル・ワース力をモデルとする引力・斥力を
受けるようにモデルされている。図８（ａ）に示すよう
に、ファン・デル・ワース力は、ｒ₀＜ｒで引力であ
り、ｒが大きくなるにつれて０に近づき、一方、０＜ｒ
＜ｒ₀で、ｒが０に近づくにつれて急速に増大するよう
な斥力を呈する。これを擬して、この発明では、ｒがｒ
₀とｒ₁の間にあるときは引力を示し、ｒがｒ₀では０
で、ｒがｒ₀以下では、次第に増加しある有限の値にと
どまるような斥力を呈するように、楕円バブル間力が定
義される。【００４８】尚、楕円バブルの方向は図８（ｃ）のよう
に一方向にそろっているとは限らず、自然現象を示すベ
クトルの方向に応じて楕円バブルの方向が異なるため、
２つの楕円バブルの中心間の距離の計算は複雑となる。
そこで、計算量を少なくするために、以下の簡略化され
た手法を用いることにより中心間距離を算出する。この
距離の算出を図９をもとに説明する。まず、図９（ａ）
のように、２つの楕円バブルの中心点を結ぶ線分を考
え、この線分の長さをｒとする。次にこの線分と楕円バ
ブルとの交点を算出し、これらの交点と楕円との中心点
の距離をそれぞれｒ_i、ｒ_jとする。ここで、この距離の
和ｒ_i＋ｒ_jの大きさを楕円バブルの中心間の距離ｒと比
較する。図９（ａ）のように、ｒ_i＋ｒ_jがｒより小さい
場合には、２つの楕円バブルは離れすぎているものと判
断し、２つの楕円バブル間に引力が発生するように定義
する。また、図９（ｂ）のように、ｒ_i＋ｒ_jがｒより大
きい場合には、２つの楕円バブルは近づきすぎているも
のと判断し、２つの楕円バブル間に斥力が発生するよう
に定義する。さらに、図９（ｃ）のように、ｒ_i＋ｒ
_{jがｒと等しい場合もしくは、実質的に等しいものと定}
_{義した所定の範囲内には、２つの楕円バブルはほぼ接し}
_{ているものと判断し、２つの楕円バブル間には力の働か}
_{ないように定義する。} _{【００４９】ｒ} ₀ は、例えば図９（ｃ）のようにして定
義されたものあり、２つの楕円バブルの間の安定な距離
を示している。さらに、ｒ₁＜ｒでは楕円バブル間力を
０と定義しているのは（図８（ｃ））、後の力学的な計
算の際に、遠隔の楕円バブルからの影響を無視して計算
量を低減するためである。【００５０】このような手法により中心間距離の計算を
簡略化することは、コンピュータの処理速度の向上及び
ベクトル場が複雑に変化する場所においても精度の高い
小さめの三角形を形成できる点で重要である。【００５１】この実施例では、楕円バブル間力ｆ（ｒ）
は、次の数１のような３次関数を使用して定義されてい
る。尚、数１では、ｒ₁＝１．５ｒ₀としている。また、
本発明はこのような楕円バブル間力の定義に限定される
ものではないが、この実施例では、数１には、楕円バブ
ルの質量が関与する項はあらわれない。従ってこの場
合、楕円バブル間力は、専ら楕円バブル間の距離と各々
の楕円バブルの直径によって決まる。【数１】【００５２】楕円バブル間力ｆ（ｒ）や楕円バブルの中
心間距離の定義は、もちろんこのような式によるものに
限定されず、この分野の当業者なら、さまざまな式を考
え出すことができるであろう。しかし、楕円バブル間力
ｆ（ｒ）は、少なくとも次のような条件は満たす必要が
ある。【００５３】(1) ｒが十分大きいとき、ｆ（ｒ）が実質
的に０になること。これは、遠方の楕円バブルからの影
響を無視できるようにするためである。 (2) ｒがある程度小さいと、引力を呈する。これは、離
れて存在する楕円バブル群を互いに近接する方向に移動
させるためである。 (3) ｒがさらに小さいと、斥力に転じる。これは、楕円
バブル群が凝集してきたとき、一定の楕円バブル間距離
を保たせるためである。実際、もし引力しか存在しない
と、楕円バブル群は、局所的に１点に収束しまい、後で
メッシュを形成することができなくなる。 (4) ｒがいくら０に近づいても、斥力は、ある一定値よ
りも大きくならない。これは、楕円バブル群の安定性を
保つためである。例えば、ファン・デル・ワース力のよ
うに、ｒが小さくなるにつれ斥力が急速に増大すると、
接近してきた楕円バブル同志が強い斥力で反発しあい、
安定な位置にとどまりにくい。【００５４】図４に戻って、先ず、ステップ３００２で
は、図２に示す頂点ノード配置手段２０２を利用して、
稜線上の頂点に楕円バブルが配置される。このとき留意
すべきなのは、非多様体データ構造においては、ここで
いう「稜線」には、面と面の間の境界としての稜線以外
にも、図７の参照番号６０１０で示すような、面上に埋
め込まれた１次元の領域も含まれることである。【００５５】次に、ステップ３００４に関連して、図２
に示す稜線ノード配置手段２０４を利用して、稜線上に
複数の楕円バブルを配置する手続きについて図１０を参
照して、説明する。【００５６】図１０において、楕円バブルを配置しよう
とする稜線は、一般的には直線とは限らず、ｓをパラメ
ータとして、Ｌ＝（ｘ（ｓ），ｙ（ｓ），ｚ（ｓ））で
あらわされるような曲線となる。そこで、例えば、パラ
メータｓが［０，１］の範囲に亙っているものとする
と、［ａ，ｂ］⊂［０，１］のような値ａ，ｂを引数と
して、上記関数Ｌを用いてそれぞれの点ａ，ｂを実際の
稜線上に投影し、その投影された点での座標（ｘ，ｙ，
ｚ）に基づき、ｄ（ｘ，ｙ，ｚ）を計算し、その値の直
径をもつ楕円バブルを、その投影された点に配置する。
こうして、ａ及びｂの２つの点にそれぞれ１つずつ楕円
バブルが配置されることになる。このような手続き（手
続き１とする）を一旦用意しておくと、先ず、０，１に
対してその手続きを適用することにより頂点への楕円バ
ブル配置が行われる。続いて、２分木的に、区間を半分
に分けてそれぞれに手続き１を再帰的に適用する手続き
２を用意しておくことにより、楕円バブルの配置が稜線
全体に亙って続けられる。【００５７】そこで、手続き１の停止条件を、配置され
ているバフルがすべて隣接するか部分的に重なっている
ことであるとすることにより、結局、与えられた稜線全
域に亙って楕円バブルが隙間なく並べられた状態が得ら
れる。【００５８】しかし、こうして単に配置された楕円バブ
ルは通常、そのまま最適な配置であることはあり得な
い。そこで、ステップ３００６以下で、楕円バブルの最
適配置を図るための動力学的シミュレーションのステッ
プに移行する。【００５９】「動力学的シミュレーション」と称する所
以は、個々の楕円バブルを、ある質量をもつ質点と見な
し（これによって、楕円バブルは、慣性モーメントをも
たないと仮定し、従って）、並進運動だけを考慮して、
上記で定義した楕円バブル間力と、粘性を考慮して、二
階の常微分方程式であるニュートン方程式を立て、この
方程式を解くことによって、楕円バブルの位置を時間的
に変化させるからである。楕円バブルがｎ個配置された
として、ニュートン方程式は、次の数２のとおりとな
る。この式からも理解されるように、楕円バブルは初期
的には、１次元の稜線上に配置されるが、運動方程式と
しては、３次元空間を想定している。このため、運動方
程式を解いていくに従って、楕円バブルが稜線上から離
れていくような動きが生じ得る。そこで、後述するよう
に、これに対処する処理がステップ３０１０で行われる
ことになる。【数２】【００６０】この数２において、ｘ_iはi番目の楕円バブ
ルのｘ座標、ｍ_iは、i番目の楕円バブルの質量、一階微
分の項は、粘性係数c_iを含む、粘性を考慮した項であ
る。１つの実施例では、質量ｍ_iは、どの楕円バブルで
も等しい値に設定される。しかし、場合によっては、質
量ｍ_iは、楕円バブルの直径または体積に比例する値に
設定される。ｙ_i、ｚ_iについても、同様である。粘性係
数c_iも、個別の楕円バブルによって異なる値をとり得る
ように表現されているが、多くの場合、それを一定値ｃ
であるとして十分である。本願発明者の試みによれば、
粘性を考慮した項を省くと、楕円バブルの振動が減衰す
ることなく残存し、安定状態に全く達することができな
い。また、右辺のｆ_xi（ｔ）、ｆ_yi（ｔ）及びｆ
_zi（ｔ）はそれぞれ、時間ｔにおける、i番目の楕円バ
ブルに対する、上記で定義した周辺の楕円バブルからの
引力・斥力のｘ，ｙ，ｚ軸成分の和である。数１から見
て取れるように、楕円バブル間力は、距離ｒ₁より遠方
では０であると定義されているので、比較的近傍の楕円
バブルからの寄与のみを考慮すればよい。また、この式
では、質量や粘性係数を個別の楕円バブルi毎に異なる
ように表現しているが、通常、個々の楕円バブルによら
ない一定値にしても差し支えない。これが、ステップ３
００６に対応する処理である。【００６１】さて、通常、多体の力学系の運動方程式は
厳密には解けないので、Runge-Kutta法などの周知の常
微分方程式の数値解析技法により、時間tをわずかな値
ずつ増分して、個々の楕円バブルの座標値を計算する。
尚、常微分方程式の数値解析技法として、本発明は、Ru
nge-Kutta法に限定されるものではなく、Adamsの方法な
どの任意の方法を使用してもよい（例えば、洲之内治
男著、サイエンス社、理工系の数学１５、「数値計
算」、１９７８年９月刊などを参照）。これが、ステッ
プ３００８に対応する処理である。【００６２】このようにして楕円バブルを移動すると、
楕円バブルは必ずしも稜線上に拘束されている訳ではな
いので、個々の楕円バブルを運動方程式に従い移動させ
てゆくに従い、楕円バブルは、稜線上から離れることが
あり得る（図１１(a)の矢印ｐ１及び楕円バブルＡ参
照）。そこで、ステップ３０１０では、１つの方法とし
て、稜線との法線方向ｐ２に従って、楕円バブルを稜線
に引き戻す、という処理を行う。運動方程式としては、
特定の楕円バブルの座標値を強制的にオフセットさせ、
そのときの楕円バブルの系の座標値と速度を初期値とし
てRunge-Kutta法の計算を再開する、という手続きをと
る。ところが、図１１の楕円バブルＢのように、矢印ｑ
１に離れたものを、法線方向である矢印ｑ２に引き戻す
と、最早それはもとの稜線上にない、という場合があり
得る。このような楕円バブルＢは破壊される。これが、
ステップ３０１０で行われる処理である。【００６３】ステップ３０１０で行われる処理として
は、次のような別の方法もある。すなわち、楕円バブル
が稜線から離れた稜線上の地点（図１１(b)ではｔ₁、こ
れを、（ｘ（ｓ₁），ｙ（ｓ₁），ｚ（ｓ₁））とする）
での接線ベクトルｑ（現在、３次元空間内での運動を考
えているので、接線ベクトルｑは３次元ベクトルであ
る）を計算し、ベクトルの移動ベクトルｐ（これは、該
稜線上の地点から、Runge-Kuttaによって計算された新
しい楕円バブルの位置までの点を結ぶものである）との
内積ｐ・ｑを計算する。さらにこれを、ｑの長さ｜ｑ｜
で割ると、それは、ほぼ、新しい楕円バブルの位置の、
ｔ₁からの、稜線に沿った距離Ｄをあらわす値となって
いる。ところが、Δｓが微小な値であるとき、以下の数
３が成立する。【数３】ΔＬ／Δｓ≒｜ｑ｜【００６４】ここで、上記Ｄは、ΔＬに相当する値と見
なすことができるので、【数４】Δｓ＝Ｄ／｜ｑ｜すなわち、 Δｓ＝ｐ・ｑ／｜ｑ｜² 【００６５】よって、このようなΔｓを以て、楕円バブ
ルを引き戻す稜線上の地点は、ｔ₂＝（ｘ（ｓ₁＋Δ
ｓ），ｙ（ｓ₁＋Δｓ），ｚ（ｓ₁＋Δｓ））と計算され
る。【００６６】次に、ステップ３０１２では、力の均衡状
態かどうかが決定される。これは、例えば次のように行
われる。すなわち、Runge-Kutta法においては（あるい
は他の常微分方程式の数値解法でも同様であるが）、各
ステップ毎に、各々の座標値の差分が計算される。そこ
で、ステップ３０１２では、各々の座標値の差分の絶対
値のうちの最大値を求め、この差分の最大値が、予め定
めた値より小さい場合に、力の均衡状態に達したと判断
する。そうでないなら、処理は、直ちにステップ３００
６に戻る。【００６７】ステップ３０１４では、個々の楕円バブル
の重なり度、すなわち楕円バブルの過疎、過密を判断す
るパラメータを算出する。すなわち、ある楕円バブルに
対して、これと隣接する全ての楕円バブルとの重なりの
距離を加算した値を算出し、この値が第１のしきい値よ
り大きい場合には、「過密」と判断して楕円バブルを減
らし、第２のしきい値より小さい場合には、「過疎」と
し判断して楕円バブルを増やす。このような制御によ
り、バブルの個数を過不足なく充填する。その算出方法
を、図１２を参照して説明する。この重なり度は、上述
の楕円バブルの中心間距離の算出方法と同様の方法を用
いて計算する。すなわち、楕円バブルＡと楕円バブルＢ
の中心点を結ぶ線分の長さをｒ、この線分と楕円バブル
Ａ，Ｂとのそれぞれの交点とそれぞれの中心までの距離
をｒ_A、ｒ_Bとすると、楕円バブルＡに対する楕円バブル
Ｂの重なり度は、｛（２ｒ_A＋ｒ_B−ｒ）／ｒ_A｝×１０
０で表される。【００６８】このような重なり度の計算が、個々の楕円
バブルについて、ステップ３０１４で行われる。従っ
て、ステップ３０１６では、個々の楕円バブルにつき、
上記重なり度が所定の閾値よりも小さく且つ上記疎ら度
が別の所定の閾値よりも大きいかについて判断が行わ
れ、もしそうなら、任意の楕円バブルの位置で適切な楕
円バブル配置が達成されている、ということであるの
で、ステップ３０２０に進み、そこで個々の楕円バブル
の中心位置に、メッシュ・ノードを配置する。もしそう
でないなら、ステップ３０１８に進み、特定の楕円バブ
ルの重なり度が上記所定の閾値よりも大きい（すなわ
ち、過密部である）という条件に応答して、その楕円バ
ブルを破壊する。また、特定の楕円バブルの疎ら度が上
記別の所定の閾値よりも小さい（すなわち、過疎部であ
る）という条件に応答して、楕円バブルを分裂させる。
実際上、楕円バブルの分裂とは、実際には、新たな楕円
バブルを、過疎であると判断された楕円バブルの近傍に
配置することを意味する。こうして楕円バブルの破壊・
分裂を行った後、ステップ３００６に戻り、Runge-Kutt
aの新たなサイクルに入る。【００６９】次に、図５を参照すると、面の頂点、稜線
及び面上に楕円バブルを配置し、動力学的モデルに従い
楕円バブルを移動させ、必要に応じて楕円バブルを破壊
・分裂させるための処理のフローチャートが示されてい
る。図５の処理は、図４の処理の完了後に行われる。こ
のとき留意すべきなのは、非多様体データ構造において
は、ここでいう「面」には、ソリッド間の境界としての
面以外にも、図７の参照番号６０１２で示すような、ソ
リッド内に埋め込まれた２次元の領域も含まれることで
ある。【００７０】図５において、ステップ４００２では、図
４のフローチャートで示したのと同様の手続きにより、
面の頂点と稜線上に、楕円バブルが最適配置される。
尚、図５のフローチャートのステップが図４のフローチ
ャートのステップの続きとして実行されるのなら、ステ
ップ４００２は既に完了しているので、スキップされる
ことになる。【００７１】次に、ステップ４００４に関連して、図２
に示す面ノード配置手段２０８を利用して、稜線上に複
数の楕円バブルを配置する手続きについて図１３を参照
して、説明する。【００７２】図１３において、楕円バブルを配置しよう
とする面は、一般的には平面とは限らず、ｕ，ｖをパラ
メータとして、ｘ＝ｘ（ｕ，ｖ）、ｙ＝ｙ（ｕ，ｖ）、
ｚ＝ｚ（ｕ，ｖ）であらわされるような曲面となる。そ
こで、例えば、パラメータｕ，ｖが［０，１］の範囲に
亙っているものとすると、パラメータｕ，ｖの［０，
１］×［０，１］の平面における矩形領域の４つの頂点
上を上記（ｘ，ｙ，ｚ）に投影した点の各々に、もしそ
の点が面領域の内部であるなら、それぞれ直径ｄ（ｘ，
ｙ，ｚ）の楕円バブルを配置するような手続きを用意し
ておき、これによって、［０，１］×［０，１］を４等
分した領域の各々にこの手続きを適用し、４等分した領
域の各々をさらに４等分してこの手続きを呼び出すとい
うことを再帰的に行い、各々の再帰呼び出しにおいて、
配置された楕円バブルの全てが互いに接触または交差し
た時点で処理をリターンするようにすることにより、与
えられた平面全体が楕円バブルで稠密に埋め尽くされる
ことになる。【００７３】尚、注意を要するのは、上記のような
［０，１］×［０，１］の矩形領域との対応づけによっ
て曲面に楕円バブルを充填するためには、その曲面は、
立方体と、位相幾何学的に同相な形状でなくてはなら
ず、例えば、図１３に示すように、その曲面の一部がト
リミングされていたり、穴があいていたりすると、楕円
バブルを配置する際に、その楕円バブルが、面領域の内
部か外部かの判定を行って、内部の場合に限って実際に
楕円バブルの配置を行う必要があることである。【００７４】以下のステップ４００６〜４０１８は、そ
れぞれ、図４とステップ３００６〜３０１８と、実質的
に同様であるが、以下のことには留意する必要がある。【００７５】先ず、ステップ４００２（あるいは図４の
フローチャートのステップ）によって最適に配置された
頂点および稜線上の（あるいは、総称的に面の境界上
の）楕円バブルは、面の内部の楕円バブルの運動方程式
の計算の際には、座標位置が不動のものとして扱われ
る。すなわち、面の内部の楕円バブルに対しては、上記
で定義された力に基づき、引力・斥力を及ぼすものとし
て運動方程式の右辺にのみあらわれ、運動方程式の左辺
の常微分の項にはあらわれない。【００７６】また、ステップ４０１０では、稜線ではな
く、面からの楕円バブルの離隔を調べる。そして、面か
ら楕円バブルが離隔していると、例えば、面の法線方向
に、楕円バブルを引き戻す。あるいは、図１１(b)で行
った方法を２次元的に適用して、ｘ方向の移動量Δｕ、
及びｙ方向の移動量Δｖをそれぞれ計算することによ
り、楕円バブルを引き戻す面上の位置を決定してもよ
い。【００７７】ステップ４０１２の力の均衡状態の判断
は、ステップ３０１２と同様の処理でよい。【００７８】ステップ４０１４での楕円バブルの重なり
度の算出は、稜線上の場合とはやや異なって、楕円バブ
ルが２次元的に隣接することになるので、図１４を参照
して説明を行う。【００７９】図１４において、楕円バブルＡを基準とし
て見たとき、その周りに６個の楕円バブルＢ、Ｃ、Ｄ、
Ｅ、Ｆ、Ｇがあり、それぞれの重なり度｛（２ｒ_a＋ｒ_b
−ｒ）／ｒ_a｝×１００＝１３０％ずつ重なりあってい
るとすると、全体の重なり度は、１３０×６＝７８０％
であると定義する。破壊するための閾値を例えば７００
％に設定してあれば、このバブルは過密と判断され、破
壊されることになる。【００８０】また、全体の重なり度が所定の下限閾値、
例えば４００％、よりも小さい場合、ステップ４０１８
で楕円バブルを新たに配置（換言すると、分裂）させる
ようにする。【００８１】こうして、ステップ４０１６では、各々の
楕円バブルの重なり度から、適切な楕円バブル配置かど
うか決定し、そうでないならステップ４０１８で、過密
部の楕円バブルを破壊し、過疎部の楕円バブルを分裂さ
せて、Runge-Kuttaの計算の次のサイクルに進み、適切
な楕円バブル配置であると判断されたなら、ステップ４
０２０で、楕円バブルの中心位置にメッシュ・ノードを
配置する。【００８２】次に、図６を参照すると、空間の境界及び
内部に楕円バブルを配置し、動力学的モデルに従い楕円
バブルを移動させ、必要に応じて楕円バブルを破壊・分
裂させるための処理のフローチャートが示されている。
図６の処理は、図４及び図５の処理の完了後に行われ
る。【００８３】図６において、ステップ５００２では、図
５のフローチャートで示したのと同様の手続きにより、
空間境界上の面の頂点と稜線上に、楕円バブルが最適配
置される。尚、図６のフローチャートのステップが図５
のフローチャートのステップの続きとして実行されるの
なら、ステップ５００２は既に完了しているので、スキ
ップされることになる。【００８４】次に、ステップ５００４に関連して、図２
に示す空間ノード配置手段２１２を利用して、空間内に
複数の楕円バブルを配置する手続きについて図１５を参
照して説明する。【００８５】図１５において、楕円バブルを配置しよう
とする空間は、一般的には立方体または直方体とは限ら
ない。そこで、立方体の８隅に楕円バブルを配置する手
続きを予め用意しておき、空間領域１７０２を十分に覆
う立方体１７０４を定義する。空間領域１７０２は、一
般的に、穴１７０６や窪みをもつことがある。さて、立
方体１７０４を８分木的に８等分し、おのおのの領域に
ついて、その候補点（すなわち、分割された立方体の８
隅）が、空間領域１７０２の内部であるかどうか判断
し、もしそうであるなら、楕円バブルをその候補点に配
置し、そうでないなら、配置しない。このような処理を
８等分された各々の立方体領域について再帰的に行う。
その再帰的処理の停止条件は、その８等分された立方体
領域において配置された全ての楕円バブルが、互いに接
触または交差することである。【００８６】以下のステップ５００６〜５０１８は、そ
れぞれ、図４とステップ３００６〜３０１８と、実質的
に同様であるが、以下のことには留意する必要がある。【００８７】先ず、ステップ５００２（あるいは図４及
び図５のフローチャートのステップ）によって最適に配
置された頂点および稜線上（総称的には、空間の境界
上）の楕円バブルは、空間の内部の楕円バブルの運動方
程式の計算の際には、座標位置が不動のものとして扱わ
れる。すなわち、空間の境界上の楕円バブルについて
は、上記で定義された力に基づき、空間内部に配置され
る楕円バブルに対して引力・斥力を及ぼすものとして運
動方程式の右辺にのみあらわれ、運動方程式の左辺の常
微分の項にはあらわれない。【００８８】また、ステップ３０１０やステップ４０１
０では、楕円バブル移動の結果として稜線や曲面から離
隔した楕円バブルを法線方向に引き戻すようにしていた
が、ステップ５０１０では、空間境界の外に移動した楕
円バブルは、単に破壊される。というのは、空間境界上
には既に、ステップ５００２（あるいは、図５のフロー
チャートの処理によって）最適に楕円バブルが配置され
ているからである。【００８９】ステップ５０１２の力の均衡状態の判断
は、ステップ３０１２またはステップ４０１２と同様の
処理でよい。【００９０】ステップ５０１４での楕円バブルの重なり
度の算出については、曲面上の場合とはやや異なって、
楕円バブルが３次元的に隣接する。この場合、最密充填
状態では、１つの楕円バブルの周囲には、１２個の楕円
バブルが隣接する。従って、３次元における楕円バブル
の重なり度及び疎ら度の閾値の適切な値の１つの例は、
それぞれ、２次元の場合の２倍の値を使用することであ
る。【００９１】こうして、ステップ５０１６では、各々の
楕円バブルの重なり度から、適切な楕円バブル配置かど
うか決定し、そうでないならステップ５０１８で、過密
部の楕円バブルを破壊し、過疎部の楕円バブルを分裂さ
せて、Runge-Kuttaの計算の次のサイクルに進み、適切
な楕円バブル配置であると判断されたなら、ステップ５
０２０で、楕円バブルの中心位置にメッシュ・ノードを
配置する。【００９２】最後に、配置されたメッシュノードからベ
クトル場により指定された方向に沿うように三角形メッ
シュを結ぶ方法について述べる。上述のように、楕円バ
ブル間力が均衡する楕円バブルの配置が求まると、こら
らの中心点を結んで三角形メッシュを生成しなけらばな
らない。従来のバブル・メッシュ法ではこの作業を計算
幾何学で研究されたデローニ三角分割法と呼ばれるアル
ゴリズムで行っていた。このアルゴリズムを用いること
により、４つの隣接する点が与えられたときに、二通り
の対角線のうちどちらを選ぶとより正三角形に近い三角
形ができるかを自動的に選択することができる。しかし
ながら、この従来のデローニ三角分割では、非等方性メ
ッシュの方向性が全く考慮されていないので本発明では
直接使うことはできない。そこで、以下に述べるよう
に、指定されたメッシュの方向性とアスペクト比を用い
て変更されたデローニ三角分割を行うことにより、指定
された方向に沿うように三角形分割を行えるようにし
た。【００９３】まず、従来のデローニ三角分割を１６図を
元に説明する。メッシュ・ノードとして与えられた４つ
の円のそれぞれの中心点に対して、２つの三角分割が考
えられるので、それぞれについて三角形の外接円を作
る。このとき、デローニ三角形、すなわち、正三角形に
より近い三角形を得るには、「外接円の中にもう一つの
点が含まれない」という条件を満たす方を選べばよい。
１６図では、メッシュ・ノードＡ、Ｂ、Ｃを結んででき
る外接円ａの中にもう一つのメッシュ・ノードＤが含ま
れないケース１がこの条件を満たしている。これが、よ
く知られているデローニ三角形の判定条件である。【００９４】これに対して、本発明では、上記の判定方
法を調べる際に、指定されている方向にアスペクト比の
逆数倍だけ縮めた座標値を用いてデローニ三角形の判定
条件を実行することにより、ベクトルにより指定された
方向に細長いメッシュができやすいように重みづけして
いる。すなわち、まず、楕円メッシュの中心点Ａ、Ｂ、
Ｃ、Ｄの結線方法を考えた場合、まずそれぞれの楕円メ
ッシュのベクトル（ベクトル場により指定されている）
を合成する。そして、この合成ベクトルを中心点の数で
割った（すなわち４で割った）大きさを有するベクトル
を加重平均ベクトルｘとして求める。そして、各中心点
の座標に加重平均ベクトルの成分の逆数倍だけ縮めた座
標値を用いて、従来と同様にデローニ三角形の判定を実
行する。そして、条件を満たす結線方法で４つの中心点
を結線することにより三角形分割ができる。例として、
２次元座標を考えた場合、加重平均ベクトルｘの方向が
（２．０、０．０）でアスペクト比が２．０と指定され
ているとすると、各中心点Ａ、Ｂ、Ｃ、Ｄのｘ座標値を
１／２．０（＝０．５）倍した座標値を用いて判定条件
を調べる。そして、デローニ三角形の判定条件を満たす
結線方法で、中心点Ａ、Ｂ、Ｃ、Ｄを結ぶ。なお、逆数
倍縮めた座標値は、あくまで判定条件を調べる場合にの
み用い、結線は本来の中心点の座標を用いている点に注
意されたい。このようにすることにより、ベクトル方向
に沿った三角形分割を行うことができる。【００９５】１７図は、従来のデローニ法と本発明の方
法による三角分割の結果を比較した図である。与えられ
たノードは正三角形のグリッドなので、従来の方法で
は、要素が正三角形になるように分割が行われている
（図１７（ａ））。これに対して、本発明の方法では、
ベクトル場の流れ方向（図１７（ｂ））が考慮された非
等方性メッシュを形成することができる（図１７
（ｃ））。【００９６】なお、以上に詳述した非等方メッシュは、
計算力学を用いた解析以外にも、コンピュータ・グラフ
ィックスやＣＡＭなどのさまざまな計算機応用分野にお
いて利用することができる。本発明の応用例として、以
下の２つを例示する。【００９７】コンピュータ・グラフィックスの応用分野
では、曲面で定義された形状を、表示のために三角形パ
ッチの集合に切り分けて近似する作業が不可欠となる。
この際、曲面を多数の細かい三角形メッシュに分割すれ
ば、このメッシュは元の曲面形状を忠実に表現できるの
で近似の質は向上するが、処理速度は低下する。逆に、
全体に粗い三角形メッシュを用いて処理速度を向上させ
ようとすると、近似精度は低下する。そこで、曲率が大
きい部分（大きく曲がっている部分）を小さい三角形で
近似し、曲率が小さい部分（平坦な部分）を大きい三角
形で近似するのが効率のよい三角形分割である。また、
曲面上の曲率の大きさは方向によって異なるので、メッ
シュの大きさを方向に応じて変えることができれば効率
を一層向上させることができる。本発明の非等方性メッ
シュ分割は、このような要求を満たすものである。例え
ば、円筒面に対しては、曲率が小さい円筒の軸方向には
長い辺を有し、曲率が大きい円周方向には短い辺を有す
る非等方性メッシュが最適な分割法となる。【００９８】また、製造過程では、ＮＣマシンによる金
型の曲面形状切削や産業用ロボットによる部品の搬送な
ど、工具や部品が三次元平面を移動して作業を行うこと
が多い。この動作の軌道（カッター・パスやロボットの
運動軌跡）を作成する際には、工具やロボットが被作業
物と干渉しないことを、コンピュータ・シミュレーショ
ンによって予め確認しておくことが必要となる。このよ
うな干渉チェックにおいては、三次元曲面の近似形状が
用いられるが、効率のよい近似形状を得るためには、本
発明の非等方性メッシュ分割を行うことが望ましい。【００９９】【発明の効果】この発明によれば、３次元図形中に１次
元の線や２次元曲面が埋めこまれて成る、いわゆる非多
様体形状データ構造に対しても、動力学的モデルに基づ
く楕円バブル充填法を適用し、最適なメッシュを自動的
に形成することが可能ならしめられる。

【図面の簡単な説明】

【図１】本発明の概要を示す図である。

【図２】本発明の処理手段のブロック図である。

【図３】本発明の全体的な処理を示すフローチャート
の図である。

【図４】稜線ノード配置手段の処理を示すフローチャ
ートの図である。

【図５】面ノード配置手段の処理を示すフローチャー
トの図である。

【図６】空間ノード配置手段の処理を示すフローチャ
ートの図である。

【図７】非多様体形状データの要素分割手順を示す図
である。

【図８】本発明で定義される楕円バブル間力を説明す
るための図である。

【図９】２つの楕円バブルの間の距離と力の関係示す
図である。

【図１０】稜線上への楕円バブル配置を示す図であ
る。

【図１１】稜線から離隔した楕円バブルを引き戻す様
子を示す図である。

【図１２】稜線上での楕円バブルの重なりを示す図で
ある。

【図１３】曲面上への楕円バブル配置を示す図であ
る。

【図１４】曲面上での楕円バブルの重なりを示す図で
ある。

【図１５】空間領域中への楕円バブル配置を示す図で
ある。

【図１６】従来のデローニ三角分割を示す図である。

【図１７】従来のデローニ法と本発明の方法による三
角分割の結果を比較した図である。

【図１８】流線のベクトル場の具体例を示す概念図で
ある。

Claims

【特許請求の範囲】

【請求項１】非等方的なメッシュを生成するメッシュ生
成方法において、所定領域内において、それぞれの要素が大きさとアスペ
クト比を有するベクトル場を指定する段階と、上記要素の大きさに対応した大きさを有し、上記アスペ
クト比に対応した方向性を有する楕円バブルを発生する
段階と、所定の安定距離よりも、上記楕円バブルが近づいた場合
には斥力を生じさせ、上記安定距離よりも、上記楕円バ
ブルが遠ざかる場合には引力を生じさせるようにバブル
間力を定義し、該バブル間力に基づいて、各々の上記楕
円バブルを移動させる段階と、上記楕円バブル同士の重なり状態を示す重なり度に応じ
て、上記楕円バブルの個数を適応的に制御する段階と、上記楕円バブル間力により上記楕円バブルを移動させ、
かつ上記楕円バブルの個数を適応的に制御することによ
り、上記楕円バブルの安定した配置を求める段階と、上記楕円バブルの中心点をもとにして、三角メッシュを
形成する段階とを有することを特徴とするメッシュ生成
方法。
【請求項２】上記ベクトル場は、上記ベクトル場は、あ
る自然現象に基づいて定められる、請求項１に記載のメ
ッシュ生成方法。
【請求項３】上記バブル間力は、２つの楕円バブルの中
心間の距離の関数である、請求項１に記載のメッシュ生
成方法。
【請求項４】上記三角メッシュを形成する段階は、上記
楕円バブルの中心点をもとにして、メッシュの方向性と
アスペクト比とを考慮したデローニ三角形の判定条件を
用いて形成されることを特徴とする請求項１乃至３に記
載のメッシュ生成方法。
【請求項５】３次元非多様体データ・モデルにおいて、
コンピュータの演算処理によってメッシュを自動生成す
る方法において、要素の大きさとアスペクト比を示すベクトル場を指定す
る段階と、上記３次元非多様体データ・モデルの全ての頂点に、コ
ンピュータによって計算される物理モデルとして、質量
を有すると共に、前記要素の大きさに対応した大きさを
有し、かつ前記アスペクト比に対応した方向を有する楕
円バブルを配置する段階と、上記３次元非多様体データ・モデルの全ての稜線上に、
コンピュータによって計算される物理モデルとして、質
量を有すると共に、前記要素の大きさに対応した大きさ
を有し、かつ前記アスペクト比に対応した方向性を有す
る複数の楕円バブルを、その外周が互いに接触または交
差するように充填する段階と、上記複数の楕円バブル間に、バブル間の距離に関連する
バブル間力を定義し、該バブル間力に基づき上記複数の
楕円バブル全体からなる力学系のニュートンの運動方程
式を解くことによって、上記稜線上の複数の楕円バブル
の移動量及び移動方向を決定し、該決定された移動量及
び移動方向に従い各々の楕円バブルを移動させる段階
と、上記各々の楕円バブルの移動量が所定の大きさよりも小
さくなったことに応答して、上記稜線上の楕円バブルの
中心点を順に結んで線分要素に分割する段階と、上記３次元非多様体データ・モデルの全ての面上の閉じ
た領域に、コンピュータによって計算される物理モデル
として、質量を有すると共に、前記要素の大きさに対応
した大きさを有し、かつ前記アスペクト比に対応した方
向を有する楕円バブルを、その外周が互いに接触または
交差するように充填する段階と、上記複数の楕円バブル間に、上記バブル間力に基づき上
記面上の複数のバブル全体からなる力学系の運動方程式
を解くことによって、上記面上の複数の楕円バブルの移
動量及び移動方向を決定し、該決定された移動量及び移
動方向に従い各々の楕円バブルを移動させる段階と、上記面上の各々の楕円バブルの移動量が所定の大きさよ
りも小さくなったことに応答して、上記面上の閉じた領
域内の楕円バブルの中心点を結んで三角形要素に分割す
る段階と、上記３次元非多様体データ・モデルの全ての空間内の閉
じた領域に、コンピュータによって計算される物理モデ
ルとして、質量を有すると共に、前記要素の大きさに対
応した大きさを有し、かつ前記アスペクト比に対応した
方向を有する楕円バブルを、その外周が互いに接触また
は交差するように充填する段階と、上記複数の楕円バブル間に、上記バブル間力に基づき上
記空間内の複数のバブル全体からなる力学系の運動方程
式を解くことによって、上記空間内の複数の楕円バブル
の移動量及び移動方向を決定し、該決定された移動量及
び移動方向に従い各々の楕円バブルを移動させる段階
と、上記空間内の各々の楕円バブルの移動量が所定の大きさ
よりも小さくなったことに応答して、上記空間内の各々
の楕円バブルの中心点を結んで四面体要素に分割する段
階とを有することを特徴とする非多様体データ・モデル
においてメッシュを自動生成する方法。
【請求項６】上記ニュートンの運動方程式には、粘性を
あらわす項が含まれてなる、請求項５に記載のメッシュ
を自動生成する方法。
【請求項７】上記ベクトル場は、ある自然現象に基づい
て定められる、請求項５に記載のメッシュを自動生成す
る方法。
【請求項８】上記バブル間力は、所定の距離よりも離隔
しているとき引力として作用し、該所定の距離よりも近
接したときは斥力として作用するように定義されてい
る、請求項５に記載のメッシュを自動生成する方法。
【請求項９】上記バブル間力は、２つの楕円バブルの中
心間の距離の関数である、請求項５に記載のメッシュを
自動生成する方法。
【請求項１０】上記三角形要素に分割する段階と、上
記四面体要素に分割する段階はそれぞれ、上記ベクトル
場の上記要素の大きさ及びアスペクト比とを考慮した２
次元と３次元のDelaunayの方法によって行われる、請求
項５に記載のメッシュを自動生成する方法。
【請求項１１】上記楕円バブル間の重なり度を決定す
る段階と、該重なり度が所定の第１の閾値よりも大きい
場合には当該バブルを破壊する段階とをさらに有する請
求項５に記載のメッシュを自動生成する方法。
【請求項１２】上記楕円バブル間の重なり度を決定す
る段階と、該重なり度が所定の第２の閾値よりも小さい
場合には、新たなバブルを発生させる段階とをさらに有
する請求項５に記載のメッシュを自動生成する方法。
【請求項１３】上記力学系の運動方程式は、Runge-Ku
ttaの方法によって数値的に解かれる、請求項５のメッ
シュを自動生成する方法。
【請求項１４】非等方的なメッシュを生成するメッシュ
生成システムにおいて、所定領域内において、ある自然現象に関係し、それぞれ
の要素が大きさとアスペクト比を有するベクトル場のデ
ータを記憶する手段と、上記要素の大きさに対応した大きさを有し、上記アスペ
クト比に対応した方向性を有する楕円バブルを表すデー
タを発生する手段と、所定の安定距離よりも、上記楕円バブルが近づいた場合
には斥力を生じさせ、上記安定距離よりも、上記楕円バ
ブルが遠ざかる場合には引力を生じさせるようにバブル
間力を定義し、該バブル間力に基づいて、各々の上記楕
円バブルデータを修正する手段と、上記楕円バブル同士の重なり状態を示す重なり度に応じ
て、上記楕円バブルの個数を適応的に制御する手段と、上記楕円バブル間力により上記楕円バブルを移動させ、
かつ上記楕円バブルの個数を適応的に制御することによ
り、上記楕円バブルの安定した配置を求める手段と、上記楕円バブルの中心点をもとにして、三角メッシュデ
ータを生成する手段とを有することを特徴とするメッシ
ュ生成システム。
【請求項１５】生成された上記三角メッシュデータに基
づいて、三角メッシュを表示する手段をさらに有するこ
とを特徴とする請求項１４に記載されたメッシュ生成シ
ステム。