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

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を参照）。残された
問題はいかにして指定されたノード間距離分布を満たし
歪みの少ない要素を生成するようなノード配置を得るか
である。特に３次元の場合は難しい。

【００１３】ところで、有限要素法などを用いた解析に
おいて、計算時間や記憶空間を増やさずに解精度を改善
するには、領域全体にわたりノード間距離（すなわち三
角形の辺の長さ）を連続的な分布関数にしたがってに変
化させなければならない。この際、解析解の値が急に変
化する領域では細かな三角要素（三角要素という用語
は、平面における三角形要素と、空間における四面体要
素の総称である）を、それ以外では粗い要素を使うこと
が大切となる。

【００１４】また、与えられたトリミング曲線上や、特
に指定された点や曲線上にノードが正確に置かれている
こと。また、解析においてこれらの曲線上で境界条件を
与えることができるように、各曲線の形状が三角形の辺
の列としてメッシュ内に表現されていなくてはならな
い。

【００１５】また、一般に、メッシュを構成する三角形
は可能な限り等辺に近いことが望ましい。なぜなら著し
く歪んだ三角形要素は解析解の精度悪化を招くからであ
る。

【００１６】さらに、大変形を扱う解析などにおいて
は、メッシュ分割すべき領域形状と要求される三角形要
素の大きさが徐々に変化することがある。この場合に
は、メッシュ分割を繰り返し行う必要があるが、この際
に新しいメッシュを全く最初から作り直すのではなく、
直前のメッシュを活用して効率良く再メッシュ分割する
ことが望まれる。また、特定の部分領域だけを局所的に
再メッシュ分割する機能も解析用モデル作成の効率化に
役立つ。しかし、上記刊行物で提示されている技法は、
必ずしもこのような要望に十分に応じ得るものではな
い。

【００１７】そこで、本願発明者は、「物理モデルによ
る自動メッシュ分割」と題する、シミュレーションとい
う雑誌の第１２巻第１号、平成５年３月、１１〜２０ペ
ージの論文において、上記の要望を実現すべく、上記ノ
ード結合法を１つの前提技術とし、理想的なメッシュ形
状を生成するような物理モデルを設計し、この物理系の
支配方程式を構築し、この支配方程式を数値的方法で解
き、その解をもとにノード配置を得、ノードをDelaunay
の方法で結んでメッシュを得る、という方策に基づき、
物理モデルとしてバブル・メッシュ分割システムを提案
した。この方法によれば、バブルは、質点をもった球状
粒子であり、その直径は、所定のノード間距離によって
与えられる。さらに、２つのバブル間には、両者が接し
ているときに最も安定な状態となるように、力の場を定
義する。また、バブルの抵抗に比例した粘性抵抗も存在
するものとする。このようなモデルを想定して、立体の
メッシュ分割を行うために、立体の頂点上に頂点バブル
を配置し、立体の稜線上に稜線バブルを配置し、立体の
面上に面バブルを配置し、空間内に空間バブルを配置
し、続いて、粘性抵抗を考慮して古典的ニュートン力学
の方程式の多体問題を解き、システム全体が所定の安定
条件を満たしたことに応答してその解を打ちきることに
よりバブルの適切な配置を決定する。この論文は、「隙
間と重なりを最小化するように最密に球（バブル）を配
置してできるパターンは、歪みが小さいDelaunay三角形
と双対なVoronoi多角形に似ている」という本願発明者
の知見に基づくものであり、本願発明も、基本的にはこ
の論文を背景技術とする。

【００１８】一方、本願発明者らは、特願昭６３−２４
４９９号明細書（特開平２−１３２５７３号公報）にお
いて、非多様体境界表現モデルにおけるソリッド・モデ
リング・システムを提示した。非多様体データ・モデル
とは、要約すると、異なる次元のソリッドの重なりや境
界を保存し、且つ表現可能とするデータ・モデルであっ
て、コンピュータ・グラフィックスにおける複数光源か
らの影の重なりを反映させる陰影表現や、建築物の構造
計算を偏微分方程式で解く際の有限要素法のためのメッ
シュを形成する場合に、形状を１次元的なソリッドと、
２次元的なソリッドと、３次元的なソリッドの和集合と
して表現する際に有利なデータ・モデルである。しか
し、上記刊行物は、非多様体データ・モデルに対して、
物理モデルによる自動メッシュ分割を適用する技法につ
いては教示するものではない。

【００１９】

【発明が解決しようとする課題】この発明の目的は、非
多様体データ・モデル上で、物理モデルによる自動メッ
シュ分割を適用するためのシステム及び方法を提供する
ことにある。

【００２０】この発明の他の目的は、非多様体データ・
モデルに複数のバブルを配置し、該複数のバブルの間で
の定義された引力・斥力に基づき力学的な方程式を解
き、バブルの中心点を結ぶことにより、自動的に最適な
三角要素分割を実行するシステム及び方法を提供するこ
とにある。

【００２１】

【課題を解決するための手段】上記目的は、次に示すよ
うな本発明のステップ、及びそれらのステップを実現す
るシステムによって達成される。

【００２２】本発明は、予め用意された非多様体データ
・モデルに対して、以下に示すステップ１乃至ステップ
４の手順により最適な数のノードを最適な位置に発生す
る。また、以下のステップ５に示したようにこのノード
発生方法をワイヤーフレーム（稜線）、サーフェス（曲
面）、ソリッド（空間領域）の順に適応し、さらに三角
要素に結んでいくことにより、３次元非多様体の自動三
角要素分割を実現するものである。

【００２３】ステップ１．ノードを直径をもったバブル
に置き換え 正三角形や正四面体に近い三角要素が得られるような理
想的なノード配置を求めるために、各ノードを直径がｄ
で与えられる球（以下バブルと呼ぶ）に置き換える。

【００２４】ステップ２．バブル間力を定義 バブル間には分子間力に似た力の場をバブル間距離の関
数として定義する。二つのバブルは接しているときが安
定状態であり、この距離よりも近づくと斥力が働き、遠
ざかると引力が働くようにする。好適なバブル間力は、
例えば三次曲線の一部分などにより定義される。

【００２５】ステップ３．動力学シミュレーションによ
るバブル間力の均衡位置の算出 要点２のバブル間力に加えて、バブルの中心点に質量と
バブルの速度に比例する粘性抵抗を与えて古典的ニュー
トン力学の運動方程式をたてる。ある初期位置から始め
てこの運動方程式をRunge-Kutta法などの数値積分によ
って解くと、時間の経過とともにバブルの配置が力の均
衡状態に近づいてゆくので、十分に近づいたところで数
値積分をうちきり最終的な解とする。ただし、ワイヤー
・フレーム上とサーフェス上のバブルはそれらの幾何要
素上からバブルが外に出ないように移動させなければな
らない。

【００２６】ステップ４．バブル個数の適応制御 要点３の動力学シミュレーションの最中にバブルどうし
の重なり具合を計算して、バブルの密度が大きすぎる場
所では適宜バブルを破壊して個数を減らし、小さすぎる
場所では分裂させて個数を増やす。このメカニズムによ
ってある領域に対して適切な数のバブルを過不足なく充
填することができる。

【００２７】ステップ５．非多様体形状上へのバブルの
配置と三角要素分割 幾何要素に対して要点１から４のようなバブルの発生方
法が用意されているとき、３次元非多様体形状の三角要
素分割は以下の手順で行うことができる。

【００２８】１．全ての頂点にバブルを配置する。 ２．全ての稜線上にバブルを最密充填する。 ３．稜線上のバブルの中心点を順に結んで線分要素に分
割する。 ４．全ての面上の閉じた領域にバブルを最密充填する。 ５．面上の閉じた領域内のバブルの中心点をDelaunayの
方法などで結んで三角形要素に分割する。 ６．全ての空間内の閉じた領域にバブルを最密充填す
る。 ７．空間内の閉じた領域内のバブルの中心点をDelaunay
の方法などで結んで四面体要素に分割する。

【００２９】尚、特に非多様体モデルにあっては、上記
バブルの充填順序は、本質的である。さらに、本発明に
あっては、非多様体モデルの通常の稜線と、次元の異な
るオブジェクト間の境界をなす稜線を区別せず同様にバ
ブル充填を行うこともまた、本質的である。

【００３０】

【実施例】以下図面に基づき、本発明の好適な実施例に
ついて説明する。

【００３１】図１は、本発明の自動三角要素分割システ
ムを実現するための処理モジュールのブロック図であ
る。この実施例では、図１に示す各々の処理モジュール
のブロックは、ＡＩＸ（ＩＢＭの商標）の下で動作する
ワークステーション上で、Ｃ言語で作成された個別のプ
ログラム、ルーチンまたは関数である。しかし、これら
個々のモジュールを、ディスクリートなハードウェアで
実現してもよいし、複数台のワークステーションをＬＡ
Ｎなどで相互接続した分散処理環境で、個別のワークス
テーションに別の処理モジュールをロードし、分散並列
処理を行うようにしてもよい。

【００３２】さて、図１において、形状入力ブロック１
は、メッシュ分割すべき入力形状を定義し、要素分割に
際して必要な情報を提供するもので、メッシュ分割すべ
き入力形状を定義する非多様体入力形状定義手段１０２
と、この形状から、これらを構成する頂点、稜線、面、
空間などの位相要素を得る入力形状位相情報抽出手段１
０４と、これらの位相要素の位置と形状などの幾何情報
を得る入力形状幾何情報抽出手段１０６とから構成され
る。

【００３３】要素分割ブロック２は、形状入力ブロック
１から順次、必要な位相的・幾何的な情報を得て、メッ
シュ・ノードを配置するとともにメッシュ要素を発生
し、さらに生成されたメッシュ出力ブロック３に順次出
力する。

【００３４】より具体的には、要素分割ブロック２は、
形状入力ブロック１からの情報に基づき頂点のノードに
バブル（詳細は後述する）を配置する頂点ノード配置手
段２０２と、形状入力ブロック１からの情報に基づき稜
線上にバブルを配置する稜線ノード配置手段２０４と、
配置された頂点及び稜線から線要素を発生する線要素発
生手段２０６と、形状入力ブロック１からの情報に基づ
き面ノード上にバブルを配置する面ノード配置手段２０
８と、配置された面ノードから２次元のDelaunayの方法
に基づき三角形要素を発生する三角形要素発生手段２１
０と、形状入力ブロック１からの情報に基づき空間ノー
ド上にバブルを配置する空間ノード発生手段２１２と、
配置された空間ノードから３次元のDelaunayの方法に基
づき四面体要素を発生する四面体要素発生手段２１４か
らなる。特に非多様体データ構造の三角形分割を適切に
行うために、要素分割ブロック２における各手段２０２
〜２１４は、図１に示す上から下に順次実行することが
本質的である。

【００３５】メッシュ出力ブロック３は、要素分割ブロ
ック２から出力されるノードと、メッシュ要素の位相的
・幾何的情報を記憶しておき、必要に応じて要素分割ブ
ロック２の処理過程に、保持していた値をフィードバッ
クする。このため、メッシュ出力ブロック３は、ノード
とメッシュ要素の位相的・幾何的情報を記憶しておくメ
ッシュ定義手段３１０と、メッシュ定義手段３１０に位
相的情報を書き込むメッシュ位相情報生成手段３０２
と、メッシュ定義手段３１０に幾何的情報を書き込むメ
ッシュ位相情報生成手段３０４と、メッシュ定義手段３
１０から情報を得て要素分割ブロック２に位相的情報を
フィードバックするメッシュ位相情報抽出手段３０６
と、メッシュ定義手段３１０から情報を得て要素分割ブ
ロック２に幾何的情報をフィードバックするメッシュ幾
何情報抽出手段３０８から構成される。

【００３６】次に、図２及び図６を参照して、本発明の
処理ステップについて説明する。図２は、本発明に係る
処理の概要ステップを示すフローチャートである。図２
では、図１の形状入力ブロック１によって、図６（ａ）
に示すような非多様体入力形状が与えられているものと
する。尚、図６（ａ）において、３次元図形である四面
体の立体６００２に、稜線６００５を境界線として２次
元図形である三角形６００４が接し、さらに三角形６０
０４に、一次元図形である直線６００８が一点６００７
を共有して接している。このような異なる次元のオブジ
ェクトからなる複合図形は、非多様体データ構造を使用
しないと表現が困難であるかまたは実質的に不可能であ
る。また、立体６００２の１つの側面にはやはり一次元
の図形である「×」形状の図形が張り付いており、さら
には、立体６００２は、自身を分ける境界線６０１２を
内在している。このような形状を保持することも、非多
様体データ構造の１つの特徴である。

【００３７】図２において、ステップ２００２では先
ず、図１の頂点ノード配置手段２０２と、稜線ノード配
置手段２０４によって、稜線上の頂点と、両端点及び内
部点と、稜線上にバブルが配置される（図６の（ｂ）及
び（ｃ）を参照）。ここで、３次元の非多様体データ構
造において注意すべきなのは、点は稜線、面及び空間の
全てにおいて境界になり得るものであり、稜線は、面及
び空間の両方において境界になり得るものであり、面
は、空間において境界になり得るものである、というこ
とである。以下の説明では、「境界」とは、この意味で
使用される。例えば、非多様体データ構造では、空間中
に、孤立した点が単独で存在し得るのであって、この場
合、その点自体が、周囲の空間との境界をなしている。

【００３８】ステップ２００４では、本願発明者による
「物理モデルによる自動メッシュ分割」と題する、シミ
ュレーションという雑誌の第１２巻第１号、平成５年３
月、１１〜２０ページの論文に述べられているような力
学的技法によってバブルが移動され、必要に応じてバブ
ルが破壊・分裂される。ステップ２００６では、バブル
位置にメッシュ・ノードが配置される。ステップ２００
７では、線要素発生手段２０６によって、バブルの中心
位置を繋ぐように線要素が発生される。尚、ステップ２
００２〜２００７は、図３のフローチャートに関連して
詳述される。

【００３９】ステップ２００８では、図１の線要素発生
手段２０６によって、線要素を発生した後（図６
（ｄ））、面ノード配置手段２０８によって面の境界上
と面上にバブルが配置され、ステップ２０１０では上記
ステップ２００４と同様の技法によってバブルが移動さ
れ、必要に応じてバブルが破壊・分裂される。ステップ
２０１２では、バブル位置にメッシュ・ノードが配置さ
れる（図６（ｅ））。すると、ステップ２０１３では、
三角形要素発生手段２１０により、２次元のDelaunayの
方法を利用して、バブルの中心位置が繋がれ、以って三
角形要素が発生される（図６（ｆ））。このとき注意す
べきであるのは、非多様体データ構造においては、図６
（ａ）の「×」形状６０１０で示されるように、２次元
の面上に１次元の線が埋めこまれていたりすることであ
る。従って、２次元のDelaunayの方法を適用する際に
は、このような埋めこまれた１次元の線に交差しないよ
うに、三角形要素を発生する必要があることに留意され
たい。ステップ２００８〜２０１３は、図４のフローチ
ャートに関連して詳述される。

【００４０】ステップ２０１４では、図１の三角形要素
発生手段２１０によって三角形要素が発生された後（図
６（ｆ））、空間ノード配置手段２１２によって空間の
境界上と空間の内部にバブルが配置され、ステップ２０
１６では上記ステップ２００４と同様の技法によってバ
ブルが移動され、必要に応じてバブルが破壊・分裂され
る。ステップ２０１８では、バブル位置にメッシュ・ノ
ードが配置され、ステップ２０１９では、３次元のDela
unayの方法を利用して、その空間内のメッシュ・ノード
に関連して、四面体要素発生手段２１４によって四面体
要素が発生される（図６（ｇ））。この場合にも、非多
様体であるが故に、空間内に埋めこまれた、空間よりも
低次元の領域に交差しないように四面体要素を発生する
必要があることに留意されたい。ステップ２０１４〜２
０１９は、図５のフローチャートに関連して詳述され
る。

【００４１】このようにして結局、図６（ｈ）に示すメ
ッシュが生成される。

【００４２】次に、図３以下を参照して、本発明に係る
メッシュ生成の詳細な処理について説明する。

【００４３】図３を参照すると、稜線上の頂点及び稜線
上にバブルを配置し、動力学的モデルに従いバブルを移
動させ、必要に応じてバブルを破壊・分裂させるための
処理のフローチャートが示されている。

【００４４】そこで先ず、本発明で使用するバブルにつ
いて説明すると、それは、直径ｄの球であって、また、
質量ｍももつような仮想的な物体である。この直径ｄ
は、一般的には、定数ではなく、ｄ≡ｄ（ｘ，ｙ，ｚ）
のように、空間の座標値に依存する関数として、メッシ
ュを最終的に利用するアプリケーションによって与えら
れる。例えば、建築物のＣＡＤデータのメッシュ生成を
行う場合、応力が大きい箇所では、比較的密なメッシュ
を生成させるように、ｄの値が小さくなる。ただし、バ
ブルは、直径ｄの剛性物体を想定している訳ではなく、
寧ろ、直径ｄというのは、バブルの近接度合い及びバブ
ル間力を示す指標として使用されるものである。従っ
て、図７（ｃ）に示すように、１つのバブルは、他のバ
ブルと、直径ｄの範囲が重なり合わないようにも存在し
得るし、互いに直径ｄの範囲が丁度接触するようにも存
在し得るし、直径ｄの範囲以内にも近接し得る。

【００４５】また、バブル同志は、分子間力としてのフ
ァン・デル・ワース力をモデルとする引力・斥力を受け
るようにモデルされている。図７（ａ）に示すように、
ファン・デル・ワース力は、ｒ₀＜ｒで引力であり、ｒ
が大きくなるにつれて０に近づき、一方、０＜ｒ＜ｒ₀
で、ｒが０に近づくにつれて急速に増大するような斥力
を呈する。これを擬して、この発明では、ｒがｒ₀とｒ₁
の間にあるときは引力を示し、ｒがｒ₀では０で、ｒが
ｒ₀以下では、次第に増加しある有限の値にとどまるよ
うな斥力を呈するように、バブル間力が定義される。
尚、このｒ₀は、例えば図８のようにして定義されたも
ので、２つのバブルの間の安定な距離を示している。さ
らに、ｒ₁＜ｒではバブル間力を０と定義しているの
は、後の力学的な計算の際に、遠隔のバブルからの影響
を無視して計算量を低減するためである。

【００４６】この実施例では、バブル間力ｆ（ｒ）は、
次の数１のような３次関数を使用して定義されている。
尚、数１では、ｒ₁＝１．５ｒ₀としている。また、本発
明はこのようなバブル間力の定義に限定されるものでは
ないが、この実施例では、数１には、バブルの質量が関
与する項はあらわれない。従ってこの場合、バブル間力
は、専らバブル間の距離と各々のバブルの直径によって
決まる。すなわち、バブルの質量は、ニュートンの運動
方程式における慣性質量としてのみ効いてくるのであっ
て、換言すれば、バブルの質量が大きい程、動きにくい
バブルである、ということになる。

【数１】

【００４７】バブル間力ｆ（ｒ）の定義は、もちろんこ
のような式によるものに限定されず、この分野の当業者
なら、さまざまな式を考え出すことができるであろう。
しかし、バブル間力ｆ（ｒ）は、少なくとも次のような
条件は満たす必要がある。

【００４８】(1) ｒが十分大きいとき、ｆ（ｒ）が実質
的に０になること。これは、遠方のバブルからの影響を
無視できるようにするためである。 (2) ｒがある程度小さいと、引力を呈する。これは、離
れて存在するバブル群を互いに近接する方向に移動させ
るためである。 (3) ｒがさらに小さいと、斥力に転じる。これは、バブ
ル群が凝集してきたとき、一定のバブル間距離を保たせ
るためである。実際、もし引力しか存在しないと、バブ
ル群は、局所的に１点に収束しまい、後でメッシュを形
成することができなくなる。 (4) ｒがいくら０に近づいても、斥力は、ある一定値よ
りも大きくならない。これは、バブル群の安定性を保つ
ためである。例えば、ファン・デル・ワース力のよう
に、ｒが小さくなるにつれ斥力が急速に増大すると、接
近してきたバブル同志が強い斥力で反発しあい、安定な
位置にとどまりにくい。

【００４９】図３に戻って、先ず、ステップ３００２で
は、図１に示す頂点ノード配置手段２０２を利用して、
稜線上の頂点にバブルが配置される。このとき留意すべ
きなのは、非多様体データ構造においては、ここでいう
「稜線」には、面と面の間の境界としての稜線以外に
も、図６の参照番号６０１０で示すような、面上に埋め
込まれた１次元の領域も含まれることである。

【００５０】次に、ステップ３００４に関連して、図１
に示す稜線ノード配置手段２０４を利用して、稜線上に
複数のバブルを配置する手続きについて図９を参照し
て、説明する。

【００５１】図９において、バブルを配置しようとする
稜線は、一般的には直線とは限らず、ｓをパラメータと
して、Ｌ＝（ｘ（ｓ），ｙ（ｓ），ｚ（ｓ））であらわ
されるような曲線となる。そこで、例えば、パラメータ
ｓが［０，１］の範囲に亙っているものとすると、
［ａ，ｂ］⊂［０，１］のような値ａ，ｂを引数とし
て、上記関数Ｌを用いてそれぞれの点ａ，ｂを実際の稜
線上に投影し、その投影された点での座標（ｘ，ｙ，
ｚ）に基づき、ｄ（ｘ，ｙ，ｚ）を計算し、その値の直
径をもつバブルを、その投影された点に配置する。こう
して、ａ及びｂの２つの点にそれぞれ１つずつバブルが
配置されることになる。このような手続き（手続き１と
する）を一旦用意しておくと、先ず、０，１に対してそ
の手続きを適用することにより頂点へのバブル配置が行
われる。続いて、２分木的に、区間を半分に分けてそれ
ぞれに手続き１を再帰的に適用する手続き２を用意して
おくことにより、バブルの配置が稜線全体に亙って続け
られる。

【００５２】そこで、手続き１の停止条件を、配置され
ているバフルがすべて隣接するか部分的に重なっている
ことであるとすることにより、結局、与えられた稜線全
域に亙ってバブルが隙間なく並べられた状態が得られ
る。

【００５３】しかし、こうして単に配置されたバブルは
通常、そのまま最適な配置であることはあり得ない。そ
こで、ステップ３００６以下で、バブルの最適配置を図
るための動力学的シミュレーションのステップに移行す
る。

【００５４】「動力学的シミュレーション」と称する所
以は、個々のバブルを、ある質量をもつ質点と見なし
（これによって、バブルは、慣性モーメントをもたない
と仮定し、従って）、並進運動だけを考慮して、上記で
定義したバブル間力と、粘性を考慮して、二階の常微分
方程式であるニュートン方程式を立て、この方程式を解
くことによって、バブルの位置を時間的に変化させるか
らである。バブルがｎ個配置されたとして、ニュートン
方程式は、次の数２のとおりとなる。この式からも理解
されるように、バブルは初期的には、１次元の稜線上に
配置されるが、運動方程式としては、３次元空間を想定
している。このため、運動方程式を解いていくに従っ
て、バブルが稜線上から離れていくような動きが生じ得
る。そこで、後述するように、これに対処する処理がス
テップ３０１０で行われることになる。

【数２】

【００５５】この数２において、ｘ_iはi番目のバブルの
ｘ座標、ｍ_iは、i番目のバブルの質量、一階微分の項
は、粘性係数c_iを含む、粘性を考慮した項である。１つ
の実施例では、質量ｍ_iは、どのバブルでも等しい値に
設定される。しかし、場合によっては、質量ｍ_iは、バ
ブルの直径または体積に比例する値に設定される。
ｙ_i、ｚ_iについても、同様である。粘性係数c_iも、個別
のバブルによって異なる値をとり得るように表現されて
いるが、多くの場合、それを一定値ｃであるとして十分
である。本願発明者の試みによれば、粘性を考慮した項
を省くと、バブルの振動が減衰することなく残存し、安
定状態に全く達することができない。また、右辺のｆ_xi
（ｔ）、ｆ_yi（ｔ）及びｆ_zi（ｔ）はそれぞれ、時間ｔ
における、i番目のバブルに対する、上記で定義した周
辺のバブルからの引力・斥力のｘ，ｙ，ｚ軸成分の和で
ある。数１から見て取れるように、バブル間力は、距離
ｒ₁より遠方では０であると定義されているので、比較
的近傍のバブルからの寄与のみを考慮すればよい。ま
た、この式では、質量や粘性係数を個別のバブルi毎に
異なるように表現しているが、通常、個々のバブルによ
らない一定値にしても差し支えない。これが、ステップ
３００６に対応する処理である。

【００５６】さて、通常、多体の力学系の運動方程式は
厳密には解けないので、Runge-Kutta法などの周知の常
微分方程式の数値解析技法により、時間tをわずかな値
ずつ増分して、個々のバブルの座標値を計算する。尚、
常微分方程式の数値解析技法として、本発明は、Runge-
Kutta法に限定されるものではなく、Adamsの方法などの
任意の方法を使用してもよい（例えば、洲之内 治男
著、サイエンス社、理工系の数学１５、「数値計算」、
１９７８年９月刊などを参照）。これが、ステップ３０
０８に対応する処理である。

【００５７】このようにしてバブルを移動すると、バブ
ルは必ずしも稜線上に拘束されている訳ではないので、
個々のバブルを運動方程式に従い移動させてゆくに従
い、バブルは、稜線上から離れることがあり得る（図１
０(a)の矢印ｐ１及びバブルＡ参照）。そこで、ステッ
プ３０１０では、１つの方法として、稜線との法線方向
ｐ２に従って、バブルを稜線に引き戻す、という処理を
行う。運動方程式としては、特定のバブルの座標値を強
制的にオフセットさせ、そのときのバブルの系の座標値
と速度を初期値としてRunge-Kutta法の計算を再開す
る、という手続きをとる。ところが、図１０のバブルＢ
のように、矢印ｑ１に離れたものを、法線方向である矢
印ｑ２に引き戻すと、最早それはもとの稜線上にない、
という場合があり得る。このようなバブルＢは破壊され
る。これが、ステップ３０１０で行われる処理である。

【００５８】ステップ３０１０で行われる処理として
は、次のような別の方法もある。すなわち、バブルが稜
線から離れた稜線上の地点（図１０(b)ではｔ₁、これ
を、（ｘ（ｓ₁），ｙ（ｓ₁），ｚ（ｓ₁））とする）で
の接線ベクトルｑ（現在、３次元空間内での運動を考え
ているので、接線ベクトルｑは３次元ベクトルである）
を計算し、ベクトルの移動ベクトルｐ（これは、該稜線
上の地点から、Runge-Kuttaによって計算された新しい
バブルの位置までの点を結ぶものである）との内積ｐ・
ｑを計算する。さらにこれを、ｑの長さ｜ｑ｜で割る
と、それは、ほぼ、新しいバブルの位置の、ｔ₁から
の、稜線に沿った距離Ｄをあらわす値となっている。と
ころが、Δｓが微小な値であるとき、以下の数３が成立
する。

【数３】ΔＬ／Δｓ≒｜ｑ｜

【００５９】ここで、上記Ｄは、ΔＬに相当する値と見
なすことができるので、

【数４】Δｓ＝Ｄ／｜ｑ｜ すなわち、 Δｓ＝ｐ・ｑ／｜ｑ｜²

【００６０】よって、このようなΔｓを以て、バブルを
引き戻す稜線上の地点は、ｔ₂＝（ｘ（ｓ₁＋Δｓ），ｙ
（ｓ₁＋Δｓ），ｚ（ｓ₁＋Δｓ））と計算される。

【００６１】次に、ステップ３０１２では、力の均衡状
態かどうかが決定される。これは、例えば次のように行
われる。すなわち、Runge-Kutta法においては（あるい
は他の常微分方程式の数値解法でも同様であるが）、各
ステップ毎に、各々の座標値の差分が計算される。そこ
で、ステップ３０１２では、各々の座標値の差分の絶対
値のうちの最大値を求め、この差分の最大値が、予め定
めた値より小さい場合に、力の均衡状態に達したと判断
する。そうでないなら、処理は、直ちにステップ３００
６に戻る。

【００６２】ステップ３０１４では、個々のバブルの重
なり度を算出する。その算出方法を、図１１を参照して
説明する。バブルＡを基準として見たとき、バブルＡと
バブルＣは、全く重なり合っていないので、バブルＡと
バブルＣの間の重なり度は、０％である。しかし、バブ
ルＢは、バブルＡの半径の１／２まで重なっているの
で、バブルＡとバブルＣの間の重なり度は、５０％であ
る。よって、バブルＡは、０＋５０＝５０％の重なり度
をもつことになり、このような指標により、バブルの重
なり度合いを検出することができる。

【００６３】一方、バブルの疎らの度合いは、図１２に
示すように、バブルＡを中心として、半径が２倍の仮想
的な球Ａ'を想定し、このＡ'内にどの程度隣接するバブ
ルが入り込んでいるかで判断する。例えば、図１２の場
合、Ａ'の幅に対してＢの半径全体が重なっているの
で、Ａ'に対するＢの重なりは、１００％であり、同様
に、Ａ'に対するＣの重なりは、７５％である。よっ
て、図１２の場合の疎ら度は、１００＋７５＝１７５％
と定義する。

【００６４】このような重なり度の計算及び疎ら度の計
算が、個々のバブルについて、ステップ３０１４で行わ
れる。従って、ステップ３０１６では、個々のバブルに
つき、上記重なり度が所定の閾値よりも小さく且つ上記
疎ら度が別の所定の閾値よりも大きいかについて判断が
行われ、もしそうなら、任意のバブルの位置で適切なバ
ブル配置が達成されている、ということであるので、ス
テップ３０２０に進み、そこで個々のバブルの中心位置
に、メッシュ・ノードを配置する。もしそうでないな
ら、ステップ３０１８に進み、特定のバブルの重なり度
が上記所定の閾値よりも大きい（すなわち、過密部であ
る）という条件に応答して、そのバブルを破壊する。過
密であると判断するための１つの目安の閾値は、１００
％である。また、特定のバブルの疎ら度が上記別の所定
の閾値よりも小さい（すなわち、過疎部である）という
条件に応答して、バブルを分裂させる。実際上、バブル
の分裂とは、実際には、新たなバブルを、過疎であると
判断されたバブルの近傍に配置することを意味する。こ
うしてバブルの破壊・分裂を行った後、ステップ３００
６に戻り、Runge-Kuttaの新たなサイクルに入る。

【００６５】次に、図４を参照すると、面の頂点、稜線
及び面上にバブルを配置し、動力学的モデルに従いバブ
ルを移動させ、必要に応じてバブルを破壊・分裂させる
ための処理のフローチャートが示されている。図４の処
理は、図３の処理の完了後に行われる。このとき留意す
べきなのは、非多様体データ構造においては、ここでい
う「面」には、ソリッド間の境界としての面以外にも、
図６の参照番号６０１２で示すような、ソリッド内に埋
め込まれた２次元の領域も含まれることである。

【００６６】図４において、ステップ４００２では、図
３のフローチャートで示したのと同様の手続きにより、
面の頂点と稜線上に、バブルが最適配置される。尚、図
４のフローチャートのステップが図３のフローチャート
のステップの続きとして実行されるのなら、ステップ４
００２は既に完了しているので、スキップされることに
なる。

【００６７】次に、ステップ４００４に関連して、図１
に示す面ノード配置手段２０８を利用して、稜線上に複
数のバブルを配置する手続きについて図１３を参照し
て、説明する。

【００６８】図１３において、バブルを配置しようとす
る面は、一般的には平面とは限らず、ｕ，ｖをパラメー
タとして、ｘ＝ｘ（ｕ，ｖ）、ｙ＝ｙ（ｕ，ｖ）、ｚ＝
ｚ（ｕ，ｖ）であらわされるような曲面となる。そこ
で、例えば、パラメータｕ，ｖが［０，１］の範囲に亙
っているものとすると、パラメータｕ，ｖの［０，１］
×［０，１］の平面における矩形領域の４つの頂点上を
上記（ｘ，ｙ，ｚ）に投影した点の各々に、もしその点
が面領域の内部であるなら、それぞれ直径ｄ（ｘ，ｙ，
ｚ）のバブルを配置するような手続きを用意しておき、
これによって、［０，１］×［０，１］を４等分した領
域の各々にこの手続きを適用し、４等分した領域の各々
をさらに４等分してこの手続きを呼び出すということを
再帰的に行い、各々の再帰呼び出しにおいて、配置され
たバブルの全てが互いに接触または交差した時点で処理
をリターンするようにすることにより、与えられた平面
全体がバブルで稠密に埋め尽くされることになる。

【００６９】尚、注意を要するのは、上記のような
［０，１］×［０，１］の矩形領域との対応づけによっ
て曲面にバブルを充填するためには、その曲面は、立方
体と、位相幾何学的に同相な形状でなくてはならず、例
えば、図１３に示すように、その曲面の一部がトリミン
グされていたり、穴があいていたりすると、バブルを配
置する際に、そのバブルが、面領域の内部か外部かの判
定を行って、内部の場合に限って実際にバブルの配置を
行う必要があることである。

【００７０】以下のステップ４００６〜４０１８は、そ
れぞれ、図３とステップ３００６〜３０１８と、実質的
に同様であるが、以下のことには留意する必要がある。

【００７１】先ず、ステップ４００２（あるいは図３の
フローチャートのステップ）によって最適に配置された
頂点および稜線上の（あるいは、総称的に面の境界上
の）バブルは、面の内部のバブルの運動方程式の計算の
際には、座標位置が不動のものとして扱われる。すなわ
ち、面の内部のバブルに対しては、上記で定義された力
に基づき、引力・斥力を及ぼすものとして運動方程式の
右辺にのみあらわれ、運動方程式の左辺の常微分の項に
はあらわれない。

【００７２】また、ステップ４０１０では、稜線ではな
く、面からのバブルの離隔を調べる。そして、面からバ
ブルが離隔していると、例えば、面の法線方向に、バブ
ルを引き戻す。あるいは、図１０(b)で行った方法を２
次元的に適用して、ｘ方向の移動量Δｕ、及びｙ方向の
移動量Δｖをそれぞれ計算することにより、バブルを引
き戻す面上の位置を決定してもよい。

【００７３】ステップ４０１２の力の均衡状態の判断
は、ステップ３０１２と同様の処理でよい。

【００７４】ステップ４０１４でのバブルの重なり度の
算出は、稜線上の場合とはやや異なって、バブルが２次
元的に隣接することになるので、図１４を参照して説明
を行う。

【００７５】図１４において、バブルＡを基準として見
たとき、その周りに６個のバブルＢ、Ｃ、Ｄ、Ｅ、Ｆ、
Ｇがあり、それぞれが半径の１／２ずつ重なりあってい
るとすると、重なり度は、５０×６＝３００％であると
定義する。この３００％は、ステップ４０１８でバブル
を破壊するための閾値の例である。

【００７６】また、図１５を参照して、疎ら度の説明を
行う。図１５において、バブルＡの周囲には、やはり６
個のバブルＢ、Ｃ、Ｄ、Ｅ、Ｆ、Ｇが、今度は重なりあ
うことなく丁度隣接して配置されていると仮定する。こ
の場合、バブルＡの直径を２倍にした仮想的な球を描い
てみると、その仮想的な球に含まれるバブルＢ、Ｃ、
Ｄ、Ｅ、Ｆ、Ｇの半径長を、バブルＡの半径で割った値
は、６００％である。そこで、この６００％を疎ら度の
下限閾値とし、それよりも疎ら度が小さい場合、ステッ
プ４０１８でバブルを新たに配置（換言すると、分裂）
させるようにする。

【００７７】こうして、ステップ４０１６では、各々の
バブルの重なり度、疎ら度から、適切なバブル配置かど
うか決定し、そうでないならステップ４０１８で、過密
部のバブルを破壊し、過疎部のバブルを分裂させて、Ru
nge-Kuttaの計算の次のサイクルに進み、適切なバブル
配置であると判断されたなら、ステップ４０２０で、バ
ブル位置にメッシュ・ノードを配置する。図１６には、
特定の平面領域において、バブルが初期配置され、Rung
e-Kuttaの計算のサイクルを続けて行くうちに、バブル
の位置（すなわち、三角要素を形成するための点の位
置）が次第に最適化されていく様子が示されている。図
１６の左側の列は、バブルの配置、右側の列は、その時
点でDelaunayの方法で形成した三角要素を示す図であ
る。図１６(a)は、バブルの初期配置、図１６(b)は、Ru
nge-Kuttaの計算サイクルが１回経過した状態の図、図
１６(c)は、Runge-Kuttaの計算サイクルが５回経過した
状態の図、図１６(d)は、Runge-Kuttaの計算サイクルが
５０回経過した状態の図である。

【００７８】次に、図５を参照すると、空間の境界及び
内部にバブルを配置し、動力学的モデルに従いバブルを
移動させ、必要に応じてバブルを破壊・分裂させるため
の処理のフローチャートが示されている。図５の処理
は、図３及び図４の処理の完了後に行われる。

【００７９】図５において、ステップ５００２では、図
４のフローチャートで示したのと同様の手続きにより、
空間境界上の面の頂点と稜線上に、バブルが最適配置さ
れる。尚、図５のフローチャートのステップが図４のフ
ローチャートのステップの続きとして実行されるのな
ら、ステップ５００２は既に完了しているので、スキッ
プされることになる。

【００８０】次に、ステップ５００４に関連して、図１
に示す空間ノード配置手段２１２を利用して、空間内に
複数のバブルを配置する手続きについて図１７を参照し
て説明する。

【００８１】図１７において、バブルを配置しようとす
る空間は、一般的には立方体または直方体とは限らな
い。そこで、立方体の８隅にバブルを配置する手続きを
予め用意しておき、空間領域１７０２を十分に覆う立方
体１７０４を定義する。空間領域１７０２は、一般的
に、穴１７０６や窪みをもつことがある。さて、立方体
１７０４を８分木的に８等分し、おのおのの領域につい
て、その候補点（すなわち、分割された立方体の８隅）
が、空間領域１７０２の内部であるかどうか判断し、も
しそうであるなら、バブルをその候補点に配置し、そう
でないなら、配置しない。このような処理を８等分され
た各々の立方体領域について再帰的に行う。その再帰的
処理の停止条件は、その８等分された立方体領域におい
て配置された全てのバブルが、互いに接触または交差す
ることである。

【００８２】以下のステップ５００６〜５０１８は、そ
れぞれ、図３とステップ３００６〜３０１８と、実質的
に同様であるが、以下のことには留意する必要がある。

【００８３】先ず、ステップ５００２（あるいは図３及
び図４のフローチャートのステップ）によって最適に配
置された頂点および稜線上（総称的には、空間の境界
上）のバブルは、空間の内部のバブルの運動方程式の計
算の際には、座標位置が不動のものとして扱われる。す
なわち、空間の境界上のバブルについては、上記で定義
された力に基づき、空間内部に配置されるバブルに対し
て引力・斥力を及ぼすものとして運動方程式の右辺にの
みあらわれ、運動方程式の左辺の常微分の項にはあらわ
れない。

【００８４】また、ステップ３０１０やステップ４０１
０では、バブル移動の結果として稜線や曲面から離隔し
たバブルを法線方向に引き戻すようにしていたが、ステ
ップ５０１０では、空間境界の外に移動したバブルは、
単に破壊される。というのは、空間境界上には既に、ス
テップ５００２（あるいは、図４のフローチャートの処
理によって）最適にバブルが配置されているからであ
る。

【００８５】ステップ５０１２の力の均衡状態の判断
は、ステップ３０１２またはステップ４０１２と同様の
処理でよい。

【００８６】ステップ５０１４でのバブルの重なり度
（あるいは、疎ら度）の算出については、曲面上の場合
とはやや異なって、バブルが３次元的に隣接する。この
場合、最密充填状態では、１つのバブルの周囲には、１
２個のバブルが隣接する。従って、３次元におけるバブ
ルの重なり度及び疎ら度の閾値の適切な値の１つの例
は、それぞれ、２次元の場合の２倍の値を使用すること
である。

【００８７】こうして、ステップ５０１６では、各々の
バブルの重なり度、疎ら度から、適切なバブル配置かど
うか決定し、そうでないならステップ５０１８で、過密
部のバブルを破壊し、過疎部のバブルを分裂させて、Ru
nge-Kuttaの計算の次のサイクルに進み、適切なバブル
配置であると判断されたなら、ステップ５０２０で、バ
ブル位置にメッシュ・ノードを配置する。

【００８８】

【発明の効果】以上説明したように、この発明によれ
ば、３次元図形中に１次元の線や２次元曲面が埋めこま
れて成る、いわゆる非多様体形状データ構造に対して
も、動力学的モデルに基づくバブル充填法を適用し、最
適なメッシュを自動的に形成することが可能ならしめら
れる。

【図面の簡単な説明】

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

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

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

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

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

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

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

【図８】 ２つのバブルの間の安定位置を示す図であ
る。

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

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

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

【図１２】 稜線上でのバブルの疎ら度の計算について
説明するための図である。

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

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

【図１５】 曲面上でのバブルの疎ら度の計算について
説明するための図である。

【図１６】 平面上で、運動方程式の計算サイクルを進
めてゆくにつれて、バブルの配置が最適化されてゆく様
子を示す図である。

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

Claims (24)

(57)【特許請求の範囲】
1. 【請求項１】３次元非多様体データ・モデルにおいて、
   コンピュータの演算処理によってメッシュを自動生成す
   る方法であって、(a) 上記３次元非多様体データ・モデ
   ルの全ての頂点に、コンピュータによって計算される物
   理モデルとして所定の大きさの直径及び質量をもつバブ
   ルを配置する段階と、(b) 上記３次元非多様体データ・
   モデルの全ての稜線上に、コンピュータによって計算さ
   れる物理モデルとして所定の大きさの直径及び質量をも
   つ複数のバブルを、その外周が互いに接触または交差す
   るように充填する段階と、(c) 上記複数のバブル間に、
   バブル間距離に関連するバブル間力を定義し、該バブル
   間力に基づき上記複数のバブル全体からなる力学系のニ
   ュートンの運動方程式を解くことによって、上記稜線上
   の複数のバブルの移動量及び移動方向を決定し、該決定
   された移動量及び移動方向に従い各々のバブルを移動さ
   せる段階と、(d) 上記各々のバブルの移動量が所定の大
   きさよりも小さくなったことに応答して、上記稜線上の
   バブルの中心点を順に結んで線分要素に分割する段階
   と、(e) 上記３次元非多様体データ・モデルの全ての面
   上の閉じた領域に、コンピュータによって計算される物
   理モデルとして所定の大きさの直径及び質量をもつバブ
   ルを、その外周が互いに接触または交差するように充填
   する段階と、(f) 上記複数のバブル間に定義されたバブ
   ル間力に基づき上記面上の複数のバブル全体からなる力
   学系の運動方程式を解くことによって、上記面上の複数
   のバブルの移動量及び移動方向を決定し、該決定された
   移動量及び移動方向に従い各々のバブルを移動させる段
   階と、(g) 上記面上の各々のバブルの移動量が所定の大
   きさよりも小さくなったことに応答して、上記面上の閉
   じた領域内のバブルの中心点を結んで三角形要素に分割
   する段階と、(h) 上記３次元非多様体データ・モデルの
   全ての空間内の閉じた領域に、コンピュータによって計
   算される物理モデルとして所定の大きさの直径及び質量
   をもつバブルを、その外周が互いに接触または交差する
   ように充填する段階と、(i) 上記複数のバブル間に定義
   されたバブル間力に基づき上記空間内の複数のバブル全
   体からなる力学系の運動方程式を解くことによって、上
   記空間内の複数のバブルの移動量及び移動方向を決定
   し、該決定された移動量及び移動方向に従い各々のバブ
   ルを移動させる段階と、(j) 上記空間内の各々のバブル
   の移動量が所定の大きさよりも小さくなったことに応答
   して、上記空間内の各々のバブルの中心点を結んで四面
   体要素に分割する段階を有する、非多様体データ・モデ
   ルにおいてメッシュを自動生成する方法。
2. 【請求項２】上記ニュートンの運動方程式には、粘性を
   あらわす項が含まれてなる、請求項１に記載のメッシュ
   を自動生成する方法。
3. 【請求項３】上記バブルは、座標位置（ｘ，ｙ，ｚ）の
   関数である直径ｄ（ｘ，ｙ，ｚ）と、（ｘ，ｙ，ｚ）に
   かかわらず一定の質量ｍをもつものである、請求項１に
   記載のメッシュを自動生成する方法。
4. 【請求項４】上記バブル間力は、予定の距離よりも離隔
   しているとき引力として作用し、該予定の距離よりも近
   接したときは斥力として作用するように定義されてい
   る、請求項１に記載のメッシュを自動生成する方法。
5. 【請求項５】上記バブル間力は、２つのバブルの直径の
   和に関連付けられた値である、請求項４に記載のメッシ
   ュを自動生成する方法。
6. 【請求項６】 上記三角形要素に分割する段階と、上記
   四面体要素に分割する段階はそれぞれ、２次元と３次元
   のDelaunayの方法によって行われる、請求項１に記載の
   メッシュを自動生成する方法。
7. 【請求項７】 上記三角形要素に分割する段階は、上記
   面内に線要素が埋めこまれている場合、該線要素に交差
   しないように上記三角形要素の線を引く段階を有する、
   請求項６に記載のメッシュを自動生成する方法。
8. 【請求項８】 上記四面体要素に分割する段階は、上記
   空間内に面要素が埋めこまれている場合、該面要素に交
   差しないように上記四面体要素の線を引く段階を有す
   る、請求項６に記載のメッシュを自動生成する方法。
9. 【請求項９】 上記稜線上の複数のバブルの移動量及び
   移動方向を決定し、該決定された移動量及び移動方向に
   従い各々のバブルを移動させる段階は、上記各々のバブ
   ルが上記稜線から離隔したかどうかを判断し、もしそう
   なら、該離隔したバブルを該稜線上に引き戻す処理をさ
   らに有する、請求項１に記載のメッシュを自動生成する
   方法。
10. 【請求項１０】 バブル間の重なり度を決定する段階
    と、該重なり度が予定の第１の閾値よりも大きいことが
    決定された場合には当該バブルを破壊する段階をさらに
    有する、請求項９に記載のメッシュを自動生成する方
    法。
11. 【請求項１１】 バブル間の隙間を決定する段階と、該
    隙間が予定の第２の閾値よりも大きいことが決定された
    場合には、該決定された位置付近に新たなバブルを発生
    させる段階をさらに有する、請求項９に記載のメッシュ
    を自動生成する方法。
12. 【請求項１２】 上記面上の複数のバブルの移動量及び
    移動方向を決定し、該決定された移動量及び移動方向に
    従い各々のバブルを移動させる段階は、上記各々のバブ
    ルが上記面から離隔したかどうかを判断し、もしそうな
    ら、該離隔したバブルを該面上に引き戻す処理をさらに
    有する、請求項１に記載のメッシュを自動生成する方
    法。
13. 【請求項１３】 バブル間の重なり度を決定する段階
    と、該重なり度が予定の第１の閾値よりも大きいことが
    決定された場合には当該バブルを破壊する段階をさらに
    有する、請求項１２に記載のメッシュを自動生成する方
    法。
14. 【請求項１４】 バブル間の隙間を決定する段階と、該
    隙間が予定の第２の閾値よりも大きいことが決定された
    場合には、該決定された位置付近に新たなバブルを発生
    させる段階をさらに有する、請求項１２に記載のメッシ
    ュを自動生成する方法。
15. 【請求項１５】 上記空間内の複数のバブルの移動量及
    び移動方向を決定し、該決定された移動量及び移動方向
    に従い各々のバブルを移動させる段階は、上記各々のバ
    ブルが上記空間内の領域から離隔したかどうかを判断
    し、もしそうなら、該離隔したバブルを破壊する段階を
    さらに有する、請求項１に記載のメッシュを自動生成す
    る方法。
16. 【請求項１６】 バブル間の重なり度を決定する段階
    と、該重なり度が予定の第１の閾値よりも大きいことが
    決定された場合には当該バブルを破壊する段階をさらに
    有する、請求項１５に記載のメッシュを自動生成する方
    法。
17. 【請求項１７】 バブル間の隙間を決定する段階と、該
    隙間が予定の第２の閾値よりも大きいことが決定された
    場合には、該決定された位置付近に新たなバブルを発生
    させる段階をさらに有する、請求項１５に記載のメッシ
    ュを自動生成する方法。
18. 【請求項１８】 上記力学系の運動方程式は、Runge-Ku
    ttaの方法によって数値的に解かれる、請求項１のメッ
    シュを自動生成する方法。
19. 【請求項１９】３次元非多様体データ・モデルにおい
    て、コンピュータの演算処理によってメッシュを自動生
    成するシステムであって、(a) 上記３次元非多様体デー
    タ・モデルの全ての頂点に、コンピュータによって計算
    される物理モデルとして所定の大きさの直径及び質量を
    もつバブルを配置する手段と、(b) 上記３次元非多様体
    データ・モデルの全ての稜線上に、コンピュータによっ
    て計算される物理モデルとして所定の大きさの直径及び
    質量をもつ複数のバブルを、その外周が互いに接触また
    は交差するように充填する手段と、(c) 上記複数のバブ
    ル間に、バブル間距離に関連して定義されたバブル間力
    に基づき上記複数のバブル全体からなる力学系のニュー
    トンの運動方程式を解くことによって、上記稜線上の複
    数のバブルの移動量及び移動方向を決定し、該決定され
    た移動量及び移動方向に従い各々のバブルを移動させる
    手段と、(d) 上記バブルの移動量が所定の大きさよりも
    小さくなったことに応答して、上記稜線上のバブルの中
    心点を順に結んで線分要素に分割する手段と、(e) 上記
    ３次元非多様体データ・モデルの全ての面上の閉じた領
    域に、コンピュータによって計算される物理モデルとし
    て所定の大きさの直径及び質量をもつバブルを、その外
    周が互いに接触または交差するように充填する手段と、
    (f) 上記複数のバブル間に定義されたバブル間力に基づ
    き上記面上の複数のバブル全体からなる力学系の運動方
    程式を解くことによって、上記面上の複数のバブルの移
    動量及び移動方向を決定し、該決定された移動量及び移
    動方向に従い各々のバブルを移動させる手段と、(g) 上
    記面上の各々のバブルの移動量が所定の大きさよりも小
    さくなったことに応答して、上記面上の閉じた領域内の
    バブルの中心点を結んで三角形要素に分割する手段と、
    (h) 上記３次元非多様体データ・モデルの全ての空間内
    の閉じた領域に、コンピュータによって計算される物理
    モデルとして所定の大きさの直径及び質量をもつバブル
    を、その外周が互いに接触または交差するように充填す
    る手段と、(i) 上記複数のバブル間に定義されたバブル
    間力に基づき上記空間内の複数のバブル全体からなる力
    学系の運動方程式を解くことによって、上記空間内の複
    数のバブルの移動量及び移動方向を決定し、該決定され
    た移動量及び移動方向に従い各々のバブルを移動させる
    手段と、(j) 上記空間内の各々のバブルの移動量が所定
    の大きさよりも小さくなったことに応答して、上記空間
    内の各々のバブルの中心点を結んで四面体要素に分割す
    る手段を有する、 非多様体データ・モデルにおいてメッシュを自動生成す
    るシステム。
20. 【請求項２０】上記ニュートンの運動方程式には、粘性
    をあらわす項が含まれてなる、請求項１９に記載のメッ
    シュを自動生成するシステム。
21. 【請求項２１】上記バブルは、座標位置（ｘ，ｙ，ｚ）
    の関数である直径ｄ（ｘ，ｙ，ｚ）をもつものである、
    請求項１９に記載のメッシュを自動生成するシステム。
22. 【請求項２２】上記バブル間力は、予定の距離よりも離
    隔しているき引力として作用し、該予定の距離よりも近
    接したときは斥力として作用するように定義されてい
    る、請求項１９に記載のメッシュを自動生成するシステ
    ム。
23. 【請求項２３】 上記三角形要素に分割する手段は、上
    記面上に存在する線要素に交差しないように２次元のDe
    launayの方法を実行する手段を有する、請求項１９に記
    載のメッシュを自動生成するシステム。
24. 【請求項２４】 上記四面体要素に分割する手段は、上
    記空間内に存在する線要素及び面要素に交差しないよう
    に３次元のDelaunayの方法を実行する手段を有する、請
    求項１９に記載のメッシュを自動生成するシステム。