JP2018067124A - シミュレーションプログラム、シミュレーション方法および情報処理装置 - Google Patents

Abstract

【課題】反復法による計算の収束を速くすることができる。【解決手段】情報処理装置１は、第１のメッシュモデル４の複数の節点の間の距離を示す距離データ７の距離に依存する閾値と第１の行列データ６ａが示す節点の組に対応する係数値とを比較して、比較結果に応じて複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新する。情報処理装置１は、第１のメッシュモデル４から、複数の評価値に基づいて選択した複数の節点のうちの一部の節点を除外した第２のメッシュモデル５ａについての第２の行列データ６ｂを生成する。そして、情報処理装置１は、第１の行列データ６ａおよび第２の行列データ６ｂを用いて、複数の節点に対応付けられた複数の変数の解を算出する。【選択図】図１

Description

本発明は、シミュレーションプログラム、シミュレーション方法および情報処理装置に関する。

構造解析や流体解析や電磁界解析などのシミュレーションにおいては、物体や空間を示す領域を複数の要素に分割したメッシュモデルを生成し、メッシュモデルから連立方程式を生成して解くことが行われている。例えば、要素の頂点である節点に対して変数を割り当てると共に、節点と節点の組に対して物理法則から導かれる係数値を割り当てる。係数値の集合である係数行列によって定義される連立方程式を解くことで、各節点における変数の値を求めることができる。連立方程式の解法として、共役勾配（ＣＧ：Conjugate Gradient）法などの反復法が用いられることがある。反復法では、近似解を求めて誤差を評価する計算を繰り返すことで、段階的に誤差を小さくして近似解を真の解に近付ける。

反復法では、代数的マルチグリッド（ＡＭＧ：Algebraic Multigrid）法を前処理として使用することがある。代数的マルチグリッド法では、オリジナルの細かいメッシュモデルから一部の節点を除外することで粗いメッシュモデルを生成する。そして、細かいメッシュモデルに対応する係数行列に加えて、粗いメッシュモデルに対応する係数行列を併用して近似解を探索する。これは、ある要素サイズをもつメッシュモデルを使って反復計算を行うと、誤差のうち当該要素サイズに応じた周波数成分が速く減衰するという性質を利用している。すなわち、粗いメッシュモデルを使用することで、細かいメッシュモデルのみでは減衰しづらい低周波成分（長波長成分）の誤差を速く減衰させることができる。

なお、例えば、構造解析や磁場解析に使用するメッシュモデルを自動生成するメッシュ生成装置が提案されている。提案のメッシュ生成装置は、物体を示す領域を複数の要素に分割し、分割した各要素のアスペクト比が所定の許容範囲に収まっているか判定する。アスペクト比が許容範囲に収まっていない長細い要素が存在する場合、メッシュ生成装置は、節点を追加することで当該長細い要素を分割する。

また、例えば、磁気デバイスの磁化分布のシミュレーションを行う磁化分布解析方法が提案されている。提案の磁化分布解析方法では、磁気デバイスのモデルを複数の要素に分割し、要素間の相対位置に基づいて係数を分類してメモリに保存する。そして、相対位置に基づいて分類された係数のうち計算に使用する係数をメモリから適宜読み出して、各要素の磁気エネルギーを合計した全磁気エネルギーの静的平衡状態を算出する。

また、例えば、有限要素法と境界要素法を用いて生成された連立方程式を解く高速演算処理方法が提案されている。提案の高速演算処理方法では、モデルを格子状の複数の要素に分割し、境界の要素が有する辺の中から同一平面上になく互いの距離が小さい２つの辺を抽出する。そして、連立方程式の係数行列の中から、抽出した２つの辺に対応する係数値を特定し、特定した係数値が係数行列の対角線近傍に位置するように並び替える。

特開平５−２９００１５号公報 特開平１０−３２５８５８号公報 国際公開第２００８／０２６２６１号

ところで、シミュレーション対象の形状やメッシュモデルの作成方法によっては、扁平形状、長細い形状、他の節点から遠い節点をもつ尖った形状などの歪んだ要素を含むメッシュモデルが作成されることがある。歪んだ要素を含むメッシュモデルを用いて代数的マルチグリッド法などを用いた反復計算を行うと、歪んだ要素を含まないメッシュモデルを用いた場合よりも計算の収束が遅くなる傾向がある。その結果、反復回数が増加する、計算時間が増加する、近似解の精度が低下するなどの問題が生じることがある。

１つの側面では、本発明は、反復法による計算の収束を速くすることができるシミュレーションプログラム、シミュレーション方法および情報処理装置を提供することを目的とする。

１つの態様では、コンピュータに以下の処理を実行させるシミュレーションプログラムが提供される。複数の節点を含む第１のメッシュモデルについて、複数の節点の間の特性を示す係数値を成分として含む第１の行列データを取得する。第１のメッシュモデルが示す複数の節点それぞれの位置座標に基づいて、複数の節点の間の距離を示す距離データを生成する。第１のメッシュモデルの中の複数通りの節点の組それぞれについて、距離データが示す節点の組に対応する距離に依存する閾値を算出し、閾値と第１の行列データが示す節点の組に対応する係数値とを比較し、比較結果に応じて複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新する。複数の評価値に基づいて複数の節点のうちの一部の節点を選択し、第１のメッシュモデルから選択した一部の節点を除外した第２のメッシュモデルについての第２の行列データを生成する。第１の行列データおよび第２の行列データを用いて、複数の節点に対応付けられた複数の変数の解を算出する。

また、１つの態様では、コンピュータが実行するシミュレーション方法が提供される。また、１つの態様では、記憶部と処理部とを有する情報処理装置が提供される。

１つの側面では、反復法による計算の収束を速くすることができる。

第１の実施の形態の情報処理装置の例を示す図である。 第２の実施の形態のシミュレーション装置のハードウェア例を示すブロック図である。 第２の実施の形態のシミュレーション装置が有する機能を表す機能ブロックを示すブロック図である。 第２の実施の形態のシミュレーション装置で実行される行列データ生成処理を示すフローチャートである。 第２の実施の形態のシミュレーション装置で実行される階層化設定処理（ｎ＝０）を示すフローチャート（その１）である。 第２の実施の形態のシミュレーション装置で実行される階層化設定処理（ｎ＝０）を示すフローチャート（その２）である。 第２の実施の形態のシミュレーション装置で実行される階層化設定処理（ｎ≧１）を示すフローチャートである。 第２の実施の形態のシミュレーション装置で実行される計算処理を示すフローチャートである。 第２の実施の形態のシミュレーション装置で実行される代数的マルチグリッド法による前処理を示すフローチャートである。 メッシュモデルの一例を示す図である。 メッシュモデルの粗視化の一例を示す図である。 磁性体メッシュモデルの一例を示す図である。 静磁ポテンシャルを求解する際の反復回数およびその計算時間を示す表である。

以下、図面を参照して実施の形態について説明する。
［第１の実施の形態］
第１の実施の形態について、図１を用いて説明する。

図１は、第１の実施の形態の情報処理装置の例を示す図である。
情報処理装置１は、解析対象とする物体や空間の領域を複数の要素（メッシュと言うことがある）に細分化して表したメッシュモデルに基づいて、構造解析や流体解析や電磁界解析などのシミュレーションを行う。情報処理装置１は、連立方程式を反復法によって解き、節点に対応付けられた変数（例えば、節点における静磁ポテンシャル）の解を求める。

情報処理装置１は、記憶部２および処理部３を有している。記憶部２は、ＲＡＭ（Random Access Memory）などの揮発性半導体メモリでもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性ストレージでもよい。処理部３は、例えば、ＣＰＵ（Central Processing Unit）やＤＳＰ（Digital Signal Processor）などのプロセッサである。ただし、処理部３は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路を含んでもよい。プロセッサは、ＲＡＭなどのメモリに記憶されたプログラムを実行する。プログラムには、後述する処理を実行させるシミュレーションプログラムが含まれる。複数のプロセッサの集合を「マルチプロセッサ」または「プロセッサ」と言うこともある。

記憶部２は、複数の節点を含む第１のメッシュモデル４について、複数の節点の間の特性を示す係数値（ａ_ij（ｉ，ｊ＝１，…，ｎ））を成分として含む第１の行列データ６ａを記憶する。第１のメッシュモデル４に含まれる要素は、格子状に細分化されていなくてもよい。第１のメッシュモデル４は、扁平形状、長細い形状、他の節点から遠い節点をもつ尖った形状などの歪んだ要素を含むことがある。第１のメッシュモデル４は、情報処理装置１によって生成されてもよいし、他の情報処理装置によって生成されてもよい。記憶部２は、第１のメッシュモデル４を記憶してもよい。第１の行列データ６ａは、例えば、第１のメッシュモデル４に含まれる節点の位置座標、適用される物理法則、解析対象とする物体の材質を示す材料データなどから生成される。第１の行列データ６ａは、情報処理装置１によって生成されてもよいし、他の情報処理装置によって作成されてもよい。

処理部３は、第１のメッシュモデル４が示す複数の節点それぞれの位置座標に基づいて、複数の節点の間の距離を示す距離データ７を生成する。距離データ７は、各節点について全ての他の節点との間の距離を網羅的に記載したものでなくてもよく、各節点について近傍にある一部の他の節点との間の距離のみを記載したものであってもよい。「距離」は、例えば、第１のメッシュモデル４上のユークリッド距離によって表現される。

処理部３は、第１のメッシュモデル４の中の複数通りの節点の組それぞれについて、距離データ７が示す節点の組に対応する距離に依存する閾値を算出する。閾値は、例えば、距離が大きいほど閾値が大きくなるように算出される。処理部３は、閾値と第１の行列データ６ａが示す節点の組に対応する係数値とを比較し、比較結果に応じて複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新する。例えば、処理部３は、係数値が閾値より大きい場合に、節点の組のうちの一方の節点に対応する評価値を増加させる。この場合、距離が小さいほど評価値は増加しやすくなり、距離が大きいほど評価値は増加しづらくなることになる。

処理部３は、複数の評価値に基づいて複数の節点のうちの一部の節点を選択し、選択した一部の節点を第１のメッシュモデル４から除外した第２のメッシュモデル５ａについての第２の行列データ６ｂを生成する。例えば、処理部３は、第１のメッシュモデル４に含まれる節点のうち、評価値が大きい節点を優先的に第２のメッシュモデル５ａに残し、評価値が大きい節点の近傍にある節点を第２のメッシュモデル５ａから除外すると判定する。そして、処理部３は、第１の行列データ６ａおよび第２の行列データ６ｂを用いて、複数の節点に対応付けられた複数の変数の解を算出する。例えば、第２の行列データ６ｂは、反復法の前処理としての代数的マルチグリッド法に使用される。

このような構成を有する情報処理装置１が実行するシミュレーション方法について説明する。
情報処理装置１では、処理部３が、まず、第１のメッシュモデル４が示す複数の節点それぞれの位置座標に基づいて、複数の節点の間の距離を示す距離データ７を生成する。例えば、メッシュモデル４において、各節点の位置座標に基づき、節点１から節点２の距離ｄ₁₂、節点２から節点３の距離ｄ₂₃など、節点ｉから節点ｊ（ｉ，ｊ（≠ｉ）＝１，…，ｎ）まで距離ｄ_ijを計算して、距離ｄ_ijを成分として含む距離データ７を生成する。

次に、処理部３は、第１のメッシュモデル４の中の複数通りの節点の組それぞれについて、距離データ７が示す節点の組に対応する距離に依存する閾値ε_ijを算出する。次に、処理部３は、閾値ε_ijと第１の行列データ６ａが示す節点の組に対応する係数値ａ_ijとを比較し、比較結果に応じて、評価データ８に含まれる複数の節点に対応する複数の評価値λ_jのうちの少なくとも一部の評価値を更新する。例えば、係数値ａ_ijが閾値ε_ijよりも大きい場合には、評価値λ_jを増加させる（例えば、λ_jを１だけ増加させる）。なお、評価データ８に含まれる各評価値の初期値は、例えば、０に設定される。

次に、処理部３は、複数の評価値λ_jに基づいて複数の節点のうちの一部の節点を選択し、選択した一部の節点を第１のメッシュモデル４から除外した第２のメッシュモデル５ａについての第２の行列データ６ｂを生成する。そして、処理部３は、第１の行列データ６ａおよび第２の行列データ６ｂを用いて、複数の節点に対応付けられた複数の変数の解を算出する。

ところで、第１のメッシュモデル４を粗視化する方法によっては、第１のメッシュモデル４から、第２のメッシュモデル５ｂが得られることがある。メッシュモデル５ｂは、節点１〜１３を含む第１のメッシュモデル４から節点３，５，９，１２が選択されて、それら以外の節点が除外されている。メッシュモデル５ｂは、メッシュモデル５ａと異なり、節点３が選択されているために、尖った形状の歪んだ要素が含まれている。このような第２のメッシュモデル５ｂに基づき、反復法を用いて連立方程式を解こうとすると、誤差の低減が遅くなり計算の収束が遅くなることがある。その結果、反復回数が増加して、計算時間が大幅に増加してしまうことがある。また、反復法によって得られる近似解の精度が低下することがある。

一方、情報処理装置１は、第１のメッシュモデル４の複数の節点の間の距離を示す距離データ７の距離に依存する閾値と第１の行列データ６ａが示す節点の組に対応する係数値とを比較して、比較結果に応じて複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新する。情報処理装置１は、第１のメッシュモデル４から、複数の評価値に基づいて選択した複数の節点のうちの一部の節点を除外した第２のメッシュモデル５ａについての第２の行列データ６ｂを生成する。そして、情報処理装置１は、第１の行列データ６ａおよび第２の行列データ６ｂを用いて、複数の節点に対応付けられた複数の変数の解を算出する。これにより、情報処理装置１は、反復法による計算の収束を速くすることができる。その結果、反復回数が減少して計算時間が減少することがある。また、反復法によって得られる近似解の精度が向上することがある。

［第２の実施の形態］
次に、第２の実施の形態を説明する。
図２は、第２の実施の形態のシミュレーション装置のハードウェア例を示すブロック図である。

シミュレーション装置１００は、ＣＰＵ１０１、ＲＡＭ１０２、ＨＤＤ１０３、画像信号処理部１０４、入力信号処理部１０５、媒体リーダ１０６および通信インタフェース１０７を有する。シミュレーション装置１００が有する上記ユニットは、バスに接続されている。なお、シミュレーション装置１００は、第１の実施の形態の情報処理装置１に対応する。ＲＡＭ１０２またはＨＤＤ１０３は、第１の実施の形態の記憶部２に対応する。ＣＰＵ１０１は、第１の実施の形態の処理部３に対応する。

ＣＰＵ１０１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ１０１は、ＨＤＤ１０３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０２にロードし、プログラムを実行する。なお、ＣＰＵ１０１は複数のプロセッサコアを備えてもよく、シミュレーション装置１００は複数のプロセッサを備えてもよく、以下で説明する処理を複数のプロセッサまたはプロセッサコアを用いて並列に実行してもよい。また、複数のプロセッサの集合（マルチプロセッサ）を「プロセッサ」と呼んでもよい。

ＲＡＭ１０２は、ＣＰＵ１０１が実行するプログラムやＣＰＵ１０１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、シミュレーション装置１００は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。

ＨＤＤ１０３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、および、データを記憶する不揮発性の記憶装置である。プログラムには、シミュレーションプログラムが含まれる。なお、シミュレーション装置１００は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。

画像信号処理部１０４は、ＣＰＵ１０１からの命令に従って、シミュレーション装置１００に接続されたディスプレイ１１１に画像を出力する。ディスプレイ１１１としては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなどの様々な種類のディスプレイを用いることができる。

入力信号処理部１０５は、シミュレーション装置１００に接続された入力デバイス１１２から入力信号を取得し、ＣＰＵ１０１に出力する。入力デバイス１１２としては、マウスやタッチパネルやタッチパッドやトラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、シミュレーション装置１００に、複数種類の入力デバイスが接続されていてもよい。

媒体リーダ１０６は、記録媒体１１３に記録されたプログラムやデータを読み取る読み取り装置である。記録媒体１１３として、例えば、磁気ディスク、光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。磁気ディスクには、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤが含まれる。光ディスクには、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）が含まれる。

媒体リーダ１０６は、例えば、記録媒体１１３から読み取ったプログラムやデータを、ＲＡＭ１０２やＨＤＤ１０３などの他の記録媒体にコピーする。読み取られたプログラムは、例えば、ＣＰＵ１０１によって実行される。なお、記録媒体１１３は可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体１１３やＨＤＤ１０３を、コンピュータ読み取り可能な記録媒体と言うことがある。

通信インタフェース１０７は、ネットワーク１１４に接続され、ネットワーク１１４を介して他の装置と通信を行うインタフェースである。通信インタフェース１０７は、スイッチなどの通信装置とケーブルで接続される有線通信インタフェースでもよいし、基地局と無線リンクで接続される無線通信インタフェースでもよい。

次に、シミュレーション装置により実現される機能を表す機能ブロックについて、図３を用いて説明する。
図３は、第２の実施の形態のシミュレーション装置が有する機能を表す機能ブロックを示すブロック図である。

シミュレーション装置１００は、物体や空間を示す領域を複数の要素に分割したメッシュモデルから式（１）で表される連立方程式を生成して解を計算する。

Ａは係数行列（メッシュ行列データ）、ｘは解ベクトル、ｂは右辺ベクトルである。
シミュレーション装置１００は、メッシュデータ保持部１１、行列データ保持部１２、ベクトルデータ保持部１３および計算結果保持部１４を有している。メッシュデータ保持部１１、行列データ保持部１２、ベクトルデータ保持部１３および計算結果保持部１４は、例えば、ＲＡＭ１０２またはＨＤＤ１０３に確保した記憶領域を用いて実装される。また、シミュレーション装置１００は、行列データ生成部１５および計算処理部１６を有している。行列データ生成部１５および計算処理部１６は、例えば、ＣＰＵ１０１が実行するプログラムモジュールを用いて実装される。

メッシュデータ保持部１１は、メッシュモデルに含まれる複数の節点および各節点の位置座標を保持する。メッシュデータ保持部１１に記憶されるメッシュデータは、シミュレーション装置１００が、ＣＡＤ（Computer Aided Design）ソフトウェアなどのメッシュ作成ソフトウェアが生成する２次元または３次元のモデルに基づいて生成してもよい。例えば、メッシュ作成ソフトウェアに対して、解析対象の表面の形状がユーザによって入力される。すると、メッシュ作成ソフトウェアが、解析対象の内部領域や周辺領域などを線によって細分化した線形式のモデルを自動生成する。この線形式のモデルには、異なる線の交点（節点）の位置座標はまだ計算されていない。シミュレーション装置１００は、シミュレーション装置１００上で生成された線形式のモデルまたは他の装置から入力された線形式のモデルを取得し、線形式のモデルから各節点の位置座標を計算して上記のメッシュデータを生成する。ただし、上記のメッシュモデルがシミュレーション装置１００の外部で生成され、シミュレーション装置１００に対して入力されてもよい。

行列データ保持部１２は、シミュレーション装置１００が生成する各種行列データを保持する。行列データは、例えば、射影行列、連立方程式の係数行列を示すメッシュ行列データ、メッシュ距離データである。

ベクトルデータ保持部１３は、シミュレーション装置１００で用いられる各種ベクトルデータを保持する。ベクトルデータは、例えば、評価ベクトルデータ、解ベクトル、右辺ベクトルである。

計算結果保持部１４は、シミュレーション装置１００で計算された結果の情報を保持する。計算結果保持部１４が保持する計算結果は、最終的な解ベクトルである。シミュレーション装置１００は、計算結果保持部１４が保持する計算結果を、ディスプレイ１１１に表示するなど様々な態様で出力デバイスに出力することができる。例えば、シミュレーション装置１００は、解ベクトルを構成する各変数の値をディスプレイ１１１に表示してもよい。また、シミュレーション装置１００は、２次元または３次元のモデル上に計算結果をマッピングして表示してもよい。静磁ポテンシャルなどの物理量を示す矢印・記号・模様などを、モデル上の各節点に対応付けて表示する（可視化する）ことが考えられる。

行列データ生成部１５は、射影行列、係数行列を生成する。また、行列データ生成部１５は、評価ベクトルデータを更新しながら、各階層で粗視化する節点の設定処理を行って、メッシュ行列データ、メッシュ距離データを生成する。０番目の階層で使用する係数行列（最も細かいメッシュモデルに対応する当初の係数行列）は、メッシュデータ保持部１１が保持するメッシュデータ中の節点数を一辺とする正方行列である。また、０番目の階層で使用する解ベクトルおよび右辺ベクトルは、メッシュデータ保持部１１が保持するメッシュデータ中の節点数を行数とする列ベクトルである。行列データ生成部１５は、メッシュデータ保持部１１が保持するメッシュデータと物理法則を示す所定の方程式とから、０番目の階層で使用する係数行列および右辺ベクトルを生成する。

その際、行列データ生成部１５は、解析対象の材質または材質に関する物性値を示す材料データを参照してもよい。材料データは、例えば、メッシュデータ保持部１１にメッシュデータと合わせて記憶される。材料データは、ユーザによって入力されてもよいし、ＣＡＤソフトウェアなどのソフトウェアによって生成されてもよい。

計算処理部１６は、シミュレーション対象の式（１）の連立方程式について、行列データ、ベクトルデータなどを用いて、計算処理を実行する。また、計算処理部１６は、このような計算処理中において、代数的マルチグリッド法による前処理を行う。

このような機能を有するシミュレーション装置１００は、まず、行列データ生成部１５が、のちの計算処理で用いられる行列データ、ベクトルデータなどを生成する。この作成の際に、行列データ生成部１５は、例えば、ｎ番目の階層を構成する節点から、粗視化により有効化されてｎ番目の階層よりも１階層下位のｎ＋１番目の階層に含まれる節点（Ｃ（Coarse）点）を設定する。さらに、行列データ生成部１５は、ｎ番目の階層を構成する節点から、粗視化によりｎ＋１番目の階層では無効化される節点（Ｆ（Fine）点）を設定する。行列データ生成部１５は各階層における行列データ、ベクトルデータなどを生成する。そして、計算処理部１６が、行列データ生成部１５により生成された行列データ、ベクトルデータを用いて、シミュレーション対象に対する計算処理を実行する。

なお、本実施の形態では、０番目の階層が最上位の階層であり、０番目から１つ下位の階層を１番目の階層とする。すなわち、ｎ番目の階層から１つ下位の階層はｎ＋１番目の階層となる。

次に、シミュレーション装置１００において、行列データ生成部１５により実行される行列データ作成処理について図４を用いて説明する。代数的マルチグリッド法を前処理として用いる反復法では、前処理を行い前処理の結果を用いて近似解を更新することが繰り返される。代数的マルチグリッド法のためのサイズの異なる複数のメッシュモデル（すなわち、複数の係数行列）は、反復計算の間に繰り返し使用される。そのため、行列データ生成部１５は、計算処理部１６が反復計算を開始する前に、以下に説明する方法で複数の係数行列を１回だけ生成すればよい。

図４は、第２の実施の形態のシミュレーション装置で実行される行列データ生成処理を示すフローチャートである。
シミュレーション装置１００の行列データ生成部１５が以下の処理を実行する。

［ステップＳ１１］ 行列データ生成部１５は、現在の階層を０番目に設定する（ｎ＝０）。
［ステップＳ１２］ 行列データ生成部１５は、メッシュ行列データとメッシュモデルの節点の間の距離とを用いて、０番目の階層に含まれる複数の節点をＣ点およびＦ点に分類する階層化設定処理（ｎ＝０）を実行する。ステップＳ１２の階層化設定処理（ｎ＝０）の詳細については後述する（図５および図６）。

［ステップＳ１３］ 行列データ生成部１５は、ｎ番目の階層のＣ点およびＦ点に基づきｎ番目の射影行列Ｐⁿを生成する。
ここで、代数的マルチグリッド法では、ある階層においてスムージングと言われる解探索処理を行って、各節点に対応させる近似値と残差を算出し、残差を１つ下の階層に伝達する。下位階層では、上位階層から伝達された残差を右辺ベクトルに設定してスムージングを行う。これを、最下層に到達するまで繰り返す。最下層では、係数行列が十分に小さいため各節点に対応させる近似値を直接法によって解き、近似値を１つ上の階層に伝達する。上位階層では、下位階層から伝達された残差に関する近似値（摂動量）を、先に算出した当該上位階層の近似値に加えた上でスムージングを行い、当該上位階層の近似値を更新する。これを、０番目の階層に到達するまで繰り返す。

上位階層から下位階層への伝達を「制限」（Restriction）と言うことがあり、下位階層から上位階層への伝達を「延長」（Prolongation）と言うことがある。上位階層の残差ベクトルを下位階層の右辺ベクトルに変換するときに、射影行列Ｐⁿの転置行列である^TＰⁿが使用される。一方、下位階層で算出された残差に関する解を上位階層にフィードバックするときに、射影行列Ｐⁿが使用される。

このような射影行列Ｐⁿは、ｎ番目の階層の節点数を行数とし、ｎ番目の階層のＣ点の節点数（すなわち、ｎ＋１番目の階層の節点数）を列数とする行列であり、以下の式（２）により作成される。

ｉはｎ番目の階層の何れかの節点を示す節点番号であり、ｊはｎ＋１番目の階層の何れかの節点を示す節点番号である。すなわち、節点ｉがＣ点の集合に含まれる場合には、Ｐ_ij＝δ_ijで表され、節点ｉがＦ点の集合に含まれる場合には、Ｐ_ij＝ｇ_ijで表される。式（２）のδ_ijは以下の式（３）で、ｇ_ijは以下の式（４）でそれぞれ表すことができる。

δ_ijは、クロネッカーのデルタである。すなわち、節点ｉと節点ｊが同じである場合には、節点ｉ，ｊの組に対応する射影行列Ｐⁿの成分は「１」になり、節点ｉと節点ｊが異なる場合には、節点ｉ，ｊの組に対応する射影行列Ｐⁿの成分は「０」になる。

式（４）のｋは、Ｆ点の集合に含まれる。
ｇ_ijは、節点ｉから節点ｊに対する依存度合に、節点ｉから他のＦ点である節点ｋを経由した節点ｊに対する依存度合を足し合わせたものを加えたものである。すなわち、粗視化により消失してしまう節点ｉによる依存度合の、節点ｊに対する重ね合わせを表している。

このような射影行列Ｐⁿを用いて上位階層から下位階層への「制限」を行うと、ある節点がＣ点である場合、その節点に対応付けられた値は下位階層でもその節点に残る。一方、ある節点がＦ点である場合、その節点に対応付けられた値は周辺のＣ点に分散して振り替えられることになる。また、このような射影行列Ｐⁿを用いて下位階層から上位階層への「延長」を行うと、下位階層の節点（Ｃ点）に対応付けられた値は上位階層でもその節点に残る。一方、下位階層に存在しない節点（Ｆ点）に対応付けるべき値は、周辺のＣ点に対応付けられた値に基づいて補完されることになる。

［ステップＳ１４］ 行列データ生成部１５は、ｎ番目の階層の係数行列Ａⁿを粗視化して、行数、列数が係数行列Ａⁿから減少した、ｎ＋１番目の階層の係数行列Ａⁿ⁺¹を生成する。

係数行列Ａⁿ⁺¹は、以下の式（５）により生成される。

なお、このようにして得られた係数行列Ａⁿ⁺¹は、ステップＳ１２で設定したＣ点の節点数を一辺とする正方行列である。
［ステップＳ１５］ 行列データ生成部１５は、係数行列Ａⁿ⁺¹の次元が、所定の最小次元であるＮ_minよりも小さいか否かを判定する。閾値としてのＮ_minは、例えば、固定値またはユーザにより与えられたものである。

次の処理は、係数行列Ａⁿ⁺¹の次元が、最小次元であるＮ_minよりも小さい場合には、行列データ作成処理を終了し、最小次元であるＮ_minよりも小さくない（Ｎ_min以上である）場合には、ステップＳ１６に進められる。

［ステップＳ１６］ 行列データ生成部１５は、現在の階層を１つ下位に下げて、ｎ＋１番目の階層に設定する。
［ステップＳ１７］ 行列データ生成部１５は、メッシュ行列データを用いて、ｎ（≧１）番目の階層に含まれる複数の節点をＣ点およびＦ点に分類する階層化設定処理（ｎ≧１）を実行する。後述するように、１番目以降の階層の階層化設定処理は、０番目の階層の階層化設定処理とは異なる。ただし、１番目以降の階層の階層化設定処理を、０番目の階層と同様の方法で行うようにしてもよい。次の処理は、再び、ステップＳ１３に進められる。ステップＳ１７の階層化設定処理（ｎ≧１）の詳細については後述する（図７および図６）。

次に、行列データ生成処理（図４）のステップＳ１２について、図５および図６を用いて説明する。
図５および図６は、第２の実施の形態のシミュレーション装置で実行される階層化設定処理（ｎ＝０）を示すフローチャートである。

シミュレーション装置１００の行列データ生成部１５が以下の処理を実行する。
［ステップＳ２１］ 行列データ生成部１５は、０番目の階層が含む節点に対応付けられた、評価ベクトルデータλの成分である評価値λ_j（ｊ＝０，１，…，Ｎ）に初期値０を設定する。ただし、Ｎは、節点の個数である。

なお、最終的な評価ベクトルデータλは、評価値λ_j（ｊ＝０，１，…，Ｎ）を成分とするｊ列の行ベクトルであり、評価値λ_jは、節点ｊに強く依存する他の節点の数を表す。

［ステップＳ２２］ 行列データ生成部１５は、０番目の階層が含む節点において、各節点の位置座標に基づき２つの節点ｉ，ｊ間の距離ｄ_ijを計算する。行列データ生成部１５は、計算した距離ｄ_ijを成分として含むメッシュ距離データを生成する。

ここで、行列データ生成部１５は、各節点について全ての他の節点との間の距離を算出しなくてもよく、各節点についてその周辺にある他の節点との間の距離を算出すればよい。例えば、行列データ生成部１５は、各節点について当該節点と同じ要素に属する他の節点との間の距離を算出する。このとき、１つの節点が複数の要素に属することがある。一例として、節点１，２，３，４が１つの四面体要素を形成し、節点１，２，５，６が１つの四面体要素を形成し、節点７，８，９，１０が１つの四面体要素を形成しているとする。節点１は２つの要素に属している。この場合、行列データ生成部１５は、節点１について、節点２，３，４，５，６との間の距離を算出することになる。節点１について、節点７，８，９，１０との間の距離は算出しなくてよい。

また、以下で使用するｉの初期値としてｉ＝０を設定する。
［ステップＳ２３］ 行列データ生成部１５は、係数行列Ａ⁰のｉ行目の非零成分の中から、ｉ＝ｊとなる対角成分を除いて、最大値の係数ｍａｘ_j≠_iａ_ijを抽出してａ^max _iとする。

［ステップＳ２４］ 行列データ生成部１５は、メッシュ距離データのｉ行目の非零成分の中から、最小距離であるｍｉｎ_iｄ_ijを抽出してｄ^min _iとする。
行列データ生成部１５は、ｊ＝０，１，…，Ｎそれぞれについて以下のステップＳ２５〜Ｓ２７を実行する。ただし、ｊはｉと異なるものとする。

［ステップＳ２５］ 行列データ生成部１５は、係数行列Ａ⁰のｉ行ｊ列の成分に対して適用するε_ijを、最小距離ｄ^min _iと、節点ｉと節点ｊの間の距離ｄ_ijとを用いて、以下の式（６）から算出する。

式（６）のεは、０＜ε＜１で設定される所定値である。εは固定値でもよいし、ユーザによって指定されてもよい。ε_ijは、距離ｄ_ijが大きいほど、大きく算出されるものである。

［ステップＳ２６］ 行列データ生成部１５は、ステップＳ２５のε_ijとステップＳ２３のａ^max _iから閾値ε_ij×ａ^max _iを算出する。行列データ生成部１５は、当該閾値と係数行列Ａ⁰のｉ行ｊ列の成分ａ_ijとを比較し、以下の式（７）を満たすか判定する。すなわち、行列データ生成部１５は、ａ_ijがε_ij×ａ^max _iより大きいか判定する。

次の処理は、式（７）を満たす場合にはステップＳ２７に進められ、式（７）を満たさない場合にはステップＳ２８に進められる。
［ステップＳ２７］ 行列データ生成部１５は、評価ベクトルデータλにおいて、式（７）を満たすｊに対応する評価値λ_jに１を加算する。

［ステップＳ２８］ 行列データ生成部１５は、ｉが節点の個数Ｎ以上か否かを判定する。
次の処理は、ｉがＮ以上である場合には、ステップＳ３０（図６）に進められ、ｉがＮ以上ではない（ｉがＮ未満である）場合には、ステップＳ２９に進められる。

［ステップＳ２９］ 行列データ生成部１５は、現在のｉに１を加算する。次の処理は、再び、ステップＳ２３に進められる。
［ステップＳ３０］ 行列データ生成部１５は、評価ベクトルデータが含む評価値λ_iのうち、まだＣ点にもＦ点にも分類されていない節点に対応する評価値の中から、最大値の評価値であるｍａｘ_iλ_iを抽出してλ^maxとする。

［ステップＳ３１］ 行列データ生成部１５は、最大の評価値λ^maxが０未満であるか否かを判定する。
次の処理は、最大の評価値λ^maxが０未満である場合には、階層化設定処理（ｎ＝０）を終了し、最大の評価値λ^maxが０未満ではない（評価値λ^maxが０以上である）場合には、ステップＳ３２に進められる。なお、まだＣ点にもＦ点にも分類されていない節点が存在しない場合には、階層化設定処理（ｎ＝０）を終了する。

［ステップＳ３２］ 行列データ生成部１５は、最大の評価値λ^maxに対応する節点の節点番号をｉ_maxとして、節点ｉ_maxをＣ点に設定する。
［ステップＳ３３］ 行列データ生成部１５は、まだＣ点にもＦ点にも分類されていない節点ｊそれぞれについて、係数行列Ａ⁰のｊ行ｉ_max列の成分ａ_jimax（ｊ＝０，１，…，Ｎ）が０であるか否か（節点ｉ_maxが節点ｊに依存するか否か）を判定する。

次の処理は、ａ_jimaxが０でないような節点ｊが存在する場合には、ステップＳ３４に進められ、そのような節点ｊが存在しない場合には、ステップＳ３６に進められる。
［ステップＳ３４］ 行列データ生成部１５は、ステップＳ３３を満たす節点ｊをＦ点に設定する。

［ステップＳ３５］ 行列データ生成部１５は、評価ベクトルデータλにおいてステップＳ３３を満たす節点ｊに対応する評価値λ_jに１を加算する。
［ステップＳ３６］ 行列データ生成部１５は、まだＣ点にもＦ点にも分類されていない節点ｊそれぞれについて、係数行列Ａ⁰のｉ_max行ｊ列の成分ａ_imaxj（ｊ＝０，１，…，Ｎ）が０であるか否か（節点ｊが節点ｉ_max（Ｃ点）に依存するか否か）を判定する。

次の処理は、ａ_imaxjが０でないような節点ｊが存在する場合には、ステップＳ３７に進められ、そのような節点ｊが存在しない場合には、ステップＳ３０に進められる。
［ステップＳ３７］ 行列データ生成部１５は、評価ベクトルデータλにおいてステップＳ３６を満たす節点ｊに対応する評価値λ_jから１を減算する。

次に、行列データ生成処理（図４）のステップＳ１７について、図７（並びに図６）を用いて説明する。
図７は、第２の実施の形態のシミュレーション装置で実行される階層化設定処理（ｎ≧１）を示すフローチャートである。

シミュレーション装置１００の行列データ生成部１５が以下の処理を実行する。
［ステップＳ４１］ 行列データ生成部１５は、ｎ番目（ｎ≧１）の階層が含む節点に対応付けられた、評価ベクトルデータλの成分である評価値λ_j（ｊ＝０，１，…，Ｎ）に初期値０を設定する。また、行列データ生成部１５は、以下で使用するｉの初期値としてｉ＝０を設定する。

［ステップＳ４２］ 行列データ生成部１５は、係数行列Ａⁿのｉ行目の非零成分の中から、ｉ＝ｊとなる対角成分を除いて、最大値の係数ｍａｘ_j≠_iａ_ijを抽出してａ^max _iとする。

行列データ生成部１５は、ｊ＝０，１，…，Ｎそれぞれについて以下のステップＳ４３，Ｓ４４を実行する。ただし、ｊはｉと異なるものとする。
［ステップＳ４３］ 行列データ生成部１５は、εとステップＳ４２のａ^max _iから閾値ε×ａ^max _iを算出する。行列データ生成部１５は、当該閾値と係数行列Ａ⁰のｉ行ｊ列の成分ａ_ijとを比較し、以下の式（８）を満たすか判定する。すなわち、行列データ生成部１５は、ａ_ijがε×ａ^max _iより大きいか判定する。

式（８）のεは、式（６）と同様に、０＜ε＜１で設定される所定値である。
次の処理は、式（８）を満たす場合には、ステップＳ４４に進められ、式（８）を満たさない場合には、ステップＳ４５に進められる。

［ステップＳ４４］ 行列データ生成部１５は、評価ベクトルデータλにおいて、式（８）を満たすｊに対応する評価値λ_jに１を加算する。
［ステップＳ４５］ 行列データ生成部１５は、ｉが節点の個数Ｎ以上か否かを判定する。

次の処理は、ｉがＮ以上である場合には、ステップＳ３０（図６）に進められ、ｉがＮ以上ではない（ｉがＮ未満である）場合には、ステップＳ４６に進められる。
なお、ステップＳ３０以降は、図６に示したフローチャートに沿った処理が実行され、各節点がＣ点またはＦ点に分類される。

［ステップＳ４６］ 行列データ生成部１５は、現在のｉに１を加算する。次の処理は、再び、ステップＳ４２に進められる。
なお、第２の実施の形態では、行列データ生成処理（図４）のステップＳ１７において、階層化設定処理（ｎ≧１）（図７および図６）を実行する場合とした。これに対し、前述のように、当該ステップＳ１７では、ｎ≧１の階層において、階層化設定処理（ｎ＝０）（図５および図６）と同様の処理を行うことも可能である。

シミュレーション装置１００では、上記のようにメッシュモデルの粗視化を行い、粗視化された階層ごとに行列データを生成し、計算処理部１６が計算処理を実行する。
次に、計算処理部１６が実行する計算処理について図８を用いて説明する。

図８は、第２の実施の形態のシミュレーション装置で実行される計算処理を示すフローチャートである。
シミュレーション装置１００の計算処理部１６が以下の処理を実行する。

［ステップＳ５１］ 計算処理部１６は、メッシュデータ保持部１１が保持するメッシュデータと物理法則を示す所定の方程式とから、式（１）の右辺ベクトルｂおよび係数行列Ａ（＝Ａ⁰）を計算する。このとき、シミュレーション対象の物体の材質を示す材料データが参照されることもある。

［ステップＳ５２］ 計算処理部１６は、解ベクトルｘの初期値ｘ₀として、０を設定する。
［ステップＳ５３］ 計算処理部１６は、以下の式（９）で表される、近似解をｘ₀とした場合の残差ｒを計算する。また、計算処理部１６は、以下で使用する変数γの初期値を０に設定し、以下で使用するベクトルｐの初期値を０に設定する。

［ステップＳ５４］ 計算処理部１６は、代数的マルチグリッド法による前処理を実行して解ベクトルｘに対応するベクトルｚを算出する。ステップＳ５４の代数的マルチグリッド法による前処理の詳細については後述する（図９）。

［ステップＳ５５］ 計算処理部１６は、変数α，β，γ，γ_oldおよびベクトルｐを式（１０）〜（１４）を用いて計算する。すなわち、計算処理部１６は、変数γを、前処理で得られたベクトルｚと直近の残差ｒとを用いて更新し、更新前の変数γの値（γ_old）と更新後の変数γの値とから変数βを更新する。そして、計算処理部１６は、前処理で得られたベクトルｚと更新した変数βとを用いてベクトルｐを更新し、ベクトルｐと係数行列Ａと変数γとを用いて変数αを更新する。

ただし、γ_old＝０の場合には、β＝０とする。

ただし、ｐ・Ａｐは、ｐとＡｐとの内積を表す。
［ステップＳ５６］ 計算処理部１６は、式（１５）に変数αおよびベクトルｐを適用して解ベクトルｘ（近似解）を更新する。

また、計算処理部１６は、式（１６）に変数α、係数行列Ａおよびベクトルｐを適用して残差ｒを更新する。

［ステップＳ５７］ 計算処理部１６は、残差ｒのノルム（絶対値）を計算する。
［ステップＳ５８］ 計算処理部１６は、所定の値である閾値δ_CGに対して、残差ｒのノルムが、δ_CG＞｜ｒ｜／｜ｂ｜であるか否か（すなわち、残差ｒが十分収束しているか否か）を判定する。

次の処理は、δ_CG＞｜ｒ｜／｜ｂ｜である（残差ｒが十分収束している）場合には、計算処理を終了し、δ_CG＞｜ｒ｜／｜ｂ｜ではない（残差ｒが十分収束していない）場合には、再び、ステップＳ５４に進められる。

次に、計算処理（図８）のステップＳ５４について、図９を用いて説明する。
図９は、第２の実施の形態のシミュレーション装置で実行される代数的マルチグリッド法による前処理を示すフローチャートである。

シミュレーション装置１００の計算処理部１６が以下の処理を実行する。
［ステップＳ６１］ 計算処理部１６は、現在の階層を０番目に設定する（ｎ＝０）。
［ステップＳ６２］ 計算処理部１６は、連立方程式Ａⁿｚⁿ＝ｒⁿについて、ｎ番目の階層の係数行列Ａⁿを利用して前進Ｇａｕｓｓ−Ｓｅｉｄｅｌ法による以下の式（１７）を少なくとも１回計算する。式（１７）を２回以上繰り返してもよい。なお、ｒ⁰（０番目の階層のｒ）はＳ５３で算出した残差ｒである。ｚの初期値は０である。

［ステップＳ６３］ 計算処理部１６は、ｎ番目の階層から１つ下位のｎ＋１番目の階層の右辺ベクトルｒⁿ⁺¹を、以下の式（１８）を用いて計算する。

［ステップＳ６４］ 計算処理部１６は、着目する階層を１つ下の階層に下げる。すなわち、計算処理部１６は、ｎに１を加算する。
［ステップＳ６５］ 計算処理部１６は、現在の階層が、所定のＮ_max番目の階層以上であるか否かを判定する。Ｎ_maxは固定値でもよいし、ユーザにより指定されてもよい。次の処理は、ｎがＮ_max以上である場合には、ステップＳ６６に進められる。また、ｎがＮ_max以上ではない（ｎがＮ_max未満である）場合には、再び、ステップＳ６２に進められる。

［ステップＳ６６］ 計算処理部１６は、連立方程式Ａⁿｚⁿ＝ｒⁿの解に相当するベクトルｚⁿを直接法で計算する。
［ステップＳ６７］ 計算処理部１６は、着目する階層を１つ上の階層に上げる。すなわち、計算処理部１６は、ｎから１を減算する。

［ステップＳ６８］ 計算処理部１６は、ｎ番目の階層の解に相当するベクトルｚⁿを、以下の式（１９）を用いて更新する。

なお、式（１９）のｎ番目の階層のベクトルｚⁿは、ベクトルｚⁿに対して、１つ下位のｎ＋１番目の階層からの寄与（ｎ＋１番目の階層のベクトルｚⁿ⁺¹を、ｎ番目の階層の射影行列Ｐⁿを用いて変換したもの）を加えることで得られる。

［ステップＳ６９］ 計算処理部１６は、連立方程式Ａⁿｚⁿ＝ｒⁿについて、ｎ番目の階層の係数行列Ａⁿを利用して後退Ｇａｕｓｓ−Ｓｅｉｄｅｌ法による以下の式（２０）を少なくとも１回計算する。式（２０）を２回以上繰り返してもよい。

［ステップＳ７０］ 計算処理部１６は、ｎが０以下であるか否かを判定する。次の処理は、ｎが０以下である場合は、代数的マルチグリッド法による前処理を終了し、ｎが０以下ではない（ｎが０より大きい）場合は、再び、ステップＳ６７に進められる。

以上により、シミュレーション装置１００では、計算処理部１６が計算処理を行って、連立方程式の解を計算する。
ここで、行列データ生成部１５がフローチャート（図５および図６）に沿って実行する、０番目の階層に含まれる節点をＣ点およびＦ点に設定する階層化設定処理（ｎ＝０）の具体例について図１０を用いて説明する。

図１０は、メッシュモデルの一例を示す図である。
メッシュモデル２０は、図１０に示されるように、４つの節点（節点１〜４）を含み、節点１〜４は０番目の階層に含まれているものとする。

このようなメッシュモデル２０は、特に、節点３が他の節点から離れて歪んだ要素を含んでいる。
まず、行列データ生成部１５は、０番目の階層が含む節点１〜４にそれぞれ対応付けられた、評価ベクトルデータλの成分である評価値λ₁〜λ₄に初期値０を設定する（ステップＳ２１）。このような評価ベクトルデータλは、以下の式（２１）のように表すことができる。

行列データ生成部１５は、０番目の階層が含む節点において、各節点の位置座標に基づき、節点ｉと節点ｊとの間の距離ｄ_ij（ｉ，ｊ＝１，２，３，４）を計算する。さらに、行列データ生成部１５は、計算した距離ｄ_ijに基づき、メッシュ距離データｄを生成する（ステップＳ２２）。このようなメッシュ距離データｄは、例えば、式（２２）のように表すことができる。

ところで、図１０の場合における係数行列Ａは、例えば、式（２３）のように表されるものとする。

行列データ生成部１５は、このような係数行列Ａの１行目の非零要素中から、対角成分を除く最大値の０．５を抽出して、ａ^max ₁とする（ステップＳ２３）。また、行列データ生成部１５は、式（２２）のメッシュ距離データｄの１行目から、対角成分を除く最小距離の１を抽出してｄ^min ₁とする（ステップＳ２４）。

行列データ生成部１５は、メッシュ距離データｄのｊ列目（ｊ＝１，２，３，４）に対する閾値ε_1jを、最小距離ｄ^min ₁（＝１）とメッシュ距離データｄとを、式（６）に適用して算出する。ただし、ここではεを０．３とした（ステップＳ２５）。

行列データ生成部１５は、ａ^max ₁に閾値ε_1jをかけると、以下の式（２４）で表される結果が得られる。

行列データ生成部１５は、式（２３）で表される係数行列Ａの１行目の対角成分以外の非零要素ａ_1j≠₁（ｊ＝２，３，４）が、式（７）をそれぞれ満たすか否かを判定する（ステップＳ２６）。

この場合、メッシュ行列データＡの１行目の（１ ０．５ ０．２ ０．１）の各列（１列除く）と、ａ^max ₁＊ε_1jの各列とをそれぞれ比較すると、メッシュ行列データＡのａ₁₂＝０．５のみがａ^max ₁＊ε_1jよりも大きい。

行列データ生成部１５は、評価ベクトルデータλのｊ＝２に対応する評価値λ₂に１を加算することで、以下の式（２５）が得られる（ステップＳ２７）。

行列データ生成部１５は、ｉ＝１は、節点の個数Ｎ＝４以上ではないことを判別し（ステップＳ２８）、メッシュ行列データＡの２〜４行目（節点２〜４）についてもそれぞれステップＳ２３〜２９の処理を繰り返し実行する。

行列データ生成部１５は、このような処理を繰り返して、以下の式（２６）で表される評価ベクトルデータλが得られる。

行列データ生成部１５は、式（２６）の評価ベクトルデータλが含む評価値から最大の評価値λ₁（＝１）を抽出して、これをλ^maxとする（ステップＳ３０）。
なお、式（２６）の評価ベクトルデータλでは、評価値λ₁，λ₂が１で最大であるものの、ここでは、評価値λ₁を抽出している。評価値λ₁に代わって、評価値λ₂を抽出しても構わない。

行列データ生成部１５は、λ^max（評価値λ₁（＝１））が０未満ではないことを判定し（ステップＳ３１）、最大の評価値λ^maxに対応する節点１をＣ点に設定する（ステップＳ３２）。

行列データ生成部１５は、係数行列Ａの非零要素のａ_j1（ｊ＝２，３，４）が０ではないことを判定し（ステップＳ３３）、節点２，３，４をＦ点に設定する（ステップＳ３３）。これにより、節点１，２，３，４がＣ点とＦ点に分類され、階層化設定処理（ｎ＝０）を終了する。

したがって、シミュレーション装置１００は、図１０のメッシュモデル２０において節点１をＣ点に、節点２〜４をＦ点に設定する。
このようにメッシュモデル２０に対して節点の間の距離を踏まえて階層化設定処理（図５および図６）を行うと、節点１がＣ点に、節点２〜４がＦ点に設定される。下位階層には、節点１がメッシュモデル２０から粗視化されて残り、メッシュモデル２０の歪んだ要素を含まない。この下位階層を用いて計算処理（図８）を行うと、反復回数の増加が抑制され、計算の収束が速くなる。

一方、メッシュモデル２０に対して、節点の間の距離を用いずに、係数行列Ａのみを用いて、歪んだ要素を除去しないような階層化設定処理（ｎ≧１）（図７および図６）によりＣ点、Ｆ点の設定を行う場合について説明する。

この場合、行列データ生成部１５が階層化設定処理（ｎ≧１）を実行すると、式（２７）で表される評価ベクトルデータλが得られる。

式（２７）の評価ベクトルデータλでは、評価値λ₃が最大であるために、節点３がＣ点として設定される。その下位階層には、メッシュモデル２０の歪みの要素である節点３がメッシュモデル２０から粗視化されて残る。この下位階層を用いて計算処理（図８）を行うと、反復回数並びに計算時間が増加してしまう。

上記シミュレーション装置１００では、メッシュ行列データと共に、節点の間の距離を用いて、メッシュモデルを粗視化するとメッシュモデルの歪みの要素を除くことができる。このような粗視化されたメッシュモデルを用いて代数的マルチグリッド法による反復計算を行うと、反復回数が減少して、計算の収束が速くなる。

また、別のメッシュモデルにおける粗視化について、図１１を用いて説明する。
図１１は、メッシュモデルの粗視化の一例を示す図である。
なお、図１１のメッシュモデル５０には、各節点の位置座標に基づき計算された、節点１から節点２の距離ｄ₁₂、節点２から節点３の距離ｄ₂₃など、節点ｉから節点ｊ（ｉ，ｊ（≠ｉ）＝１，…，ｎ）まで距離ｄ_ijが付されている。

図１１に示すメッシュモデル５０に対して、上記で説明した階層化設定処理（ｎ＝０）（図５および図６）を実行して、メッシュモデル５０に含まれる複数の節点からＣ点およびＦ点を設定する。

このようにしてメッシュモデル５０から粗視化されたメッシュモデル６０ａが得られる。メッシュモデル６０ａは、メッシュモデル５０の節点から、節点５，９，１２がＣ点として設定されて、それら以外の節点はＦ点として設定されている。このようなメッシュモデル５０ａに含まれる節点５，９，１２には、細長い形状、尖った形状などの歪んだ要素が含まれていない。したがって、メッシュモデル６０ａを用いて代数的マルチグリッド法による反復計算を行うと、反復回数が減少して、計算速度が速くなる。

一方、メッシュモデル５０を粗視化する方法によっては、メッシュモデル５０から、メッシュモデル６０ｂが得られる。メッシュモデル６０ｂは、メッシュモデル５０から節点３，５，９，１２が選択されて、それら以外の節点が除外されている。メッシュモデル６０ｂは、メッシュモデル６０ａと異なり、節点３が選択されているために、尖った形状の歪んだ要素が含まれている。このようなメッシュモデル６０ｂに基づき、反復法を用いて、複数の変数の解を算出しようとすると、反復回数が増加して、計算時間が大幅に増加してしまう。

ここで、上記を踏まえて、シミュレーション装置１００が実行するシミュレーションにより静磁界計算を行う場合について説明する。
静磁界とは、磁性体自らの時間的に定常な磁化により生成される磁界であり、その値は、静磁ポテンシャルと呼ばれる磁界のエネルギーの空間微分として与えられる。静磁ポテンシャルは、マクスウェル方程式から導出される式（２８）で求められる。

なお、φは、静磁ポテンシャル、Ｍは、磁性体の平均磁化である。
数値計算上は静磁ポテンシャルを求める際、微分方程式を離散化してメッシュモデルの節点を自由度とする連立方程式の解として計算を行う。この際、式（２８）の式（Ａ）との対応関係から、式（２８）の左辺から係数行列Ａが計算され、式（２８）の右辺から右辺ベクトルｂが計算される。

この場合の磁性体メッシュモデルについて図１２を用いて説明する。
図１２は、磁性体メッシュモデルの一例を示す図である。
なお、図１２では、磁性体メッシュモデル３０に対してｘ−ｚ平面であって、ｘ＞０およびｚ＞０の領域について示されている。

磁性体メッシュモデル３０では、立方体状であって、例えば、ネオジム磁石といった磁性体３１が空気領域３２の中心部（原点）に配置されており、磁性体３１から空気領域３２に静磁場が生じている。

磁性体メッシュモデル３０では、空気領域３２の磁性体３１に接する面のメッシュ幅（例えば、幅ｗ０）を細かく設定し、さらに、外部境界側の面のメッシュ幅（例えば、幅ｗ１）を非常に粗く設定している。このため、空気領域３２では、メッシュ内の最小の辺の長さに対する最大の辺の長さの比が１０を超える歪んだ形状になっている。

次に、磁性体メッシュモデル３０に関して式（２８）の連立方程式を解いて計算結果が１回得られるまでに要する反復回数および計算時間について、図１３を用いて説明する。
図１３は、静磁ポテンシャルを求解する際の反復回数およびその計算時間を示す表である。

磁性体メッシュモデル３０の要素数を８０，０００〜５，２００，０００の間に設定し、式（２８）の連立方程式を解くまでに要する反復回数および計測時間の計測を行った。
表４０は、当該計測結果について、式（２８）の連立方程式を解くための磁性体メッシュモデル３０の粗視化の際に、メッシュの歪み要素を除去していない場合（歪み要素非除去）を示している。また、表４０は、当該計測結果について、式（２８）の連立方程式を解くための磁性体メッシュモデル３０の粗視化の際に、メッシュ行列データと共に節点間の距離を利用して、メッシュの歪み要素を除去した場合（歪み要素除去）を示している。

この表４０によれば、各行列の次元（要素数）において、反復回数は、歪み要素除去の場合の方が、歪み要素非除去の場合に比べて、約３５％から約５０％減少している。これに伴い、計算時間も同様に、約３３％から約５０％減少している。特に、行列の次元が大きいほど、反復回数および計算時間がより改善されている。

上記の処理機能は、コンピュータによって実現することができる。その場合、コンピュータが有すべき機能の処理内容を記述したプログラムが提供される。そのプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。

処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体には、磁気記憶装置、光ディスク、光磁気記憶媒体、半導体メモリなどがある。磁気記憶装置には、ＨＤＤ、フレキシブルディスク（ＦＤ）、磁気テープ（ＭＴ）などがある。光ディスクには、ＤＶＤ、ＤＶＤ−ＲＡＭ、ＣＤ−ＲＯＭ（Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）などがある。光磁気記憶媒体には、ＭＯなどがある。半導体メモリには、ＵＳＢ（Universal Serial Bus）メモリなどのフラッシュメモリがある。

上記プログラムを流通させる場合には、例えば、そのプログラムが記録されたＤＶＤ、ＣＤ−ＲＯＭなどの可搬型記録媒体が販売される。また、プログラムをサーバコンピュータに格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することもできる。

上記プログラムを実行するコンピュータは、例えば、可搬型記録媒体に記録されたプログラム若しくはサーバコンピュータから転送されたプログラムを、自己の記憶装置に格納する。そして、コンピュータは、自己の記憶装置からプログラムを読み取り、プログラムに従った処理を実行する。なお、コンピュータは、可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することもできる。また、コンピュータは、サーバコンピュータからプログラムが転送されるごとに、逐次、受け取ったプログラムに従った処理を実行することもできる。

また、プログラムで記述された処理の一部または全てを、電子回路に置き換えることが可能である。例えば、上記の処理機能の少なくとも一部を、ＤＳＰ、ＡＳＩＣ、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。

なお、上述の実施の形態は、実施の形態の要旨を逸脱しない範囲内において種々の変更を加えることができる。
さらに、上述の実施の形態は、多数の変形、変更が当業者にとって可能であり、説明した正確な構成及び応用例に限定されるものではない。

１ 情報処理装置
２ 記憶部
３ 処理部
４ 第１のメッシュモデル
５ａ，５ｂ 第２のメッシュモデル
６ａ 第１の行列データ
６ｂ 第２の行列データ
７ 距離データ
８ 評価データ

Claims (5)

1. コンピュータに、
   複数の節点を含む第１のメッシュモデルについて、前記複数の節点の間の特性を示す係数値を成分として含む第１の行列データを取得し、
   前記第１のメッシュモデルが示す前記複数の節点それぞれの位置座標に基づいて、前記複数の節点の間の距離を示す距離データを生成し、
   前記第１のメッシュモデルの中の複数通りの節点の組それぞれについて、前記距離データが示す前記節点の組に対応する距離に依存する閾値を算出し、前記閾値と前記第１の行列データが示す前記節点の組に対応する係数値とを比較し、比較結果に応じて前記複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新し、
   前記複数の評価値に基づいて前記複数の節点のうちの一部の節点を選択し、前記第１のメッシュモデルから前記選択した一部の節点を除外した第２のメッシュモデルについての第２の行列データを生成し、
   前記第１の行列データおよび前記第２の行列データを用いて、前記複数の節点に対応付けられた複数の変数の解を算出する、
   処理を実行させるシミュレーションプログラム。
2. 前記節点の組に対応する距離が大きいほど前記閾値が大きく算出され、
   前記少なくとも一部の評価値の更新では、前記節点の組に対応する係数値が前記閾値より大きい場合に、前記節点の組の一方の節点に対応する評価値を増加させる、
   請求項１記載のシミュレーションプログラム。
3. 前記節点の組は第１の節点と第２の節点との組であり、
   前記閾値の算出では、前記第２の節点を含む２以上の節点それぞれと前記第１の節点との間の距離の中から最小距離を特定し、前記第１の節点と前記第２の節点との間の距離および前記最小距離に基づいて前記閾値を算出する、
   請求項１記載のシミュレーションプログラム。
4. コンピュータが実行するシミュレーション方法であって、
   複数の節点を含む第１のメッシュモデルについて、前記複数の節点の間の特性を示す係数値を成分として含む第１の行列データを取得し、
   前記第１のメッシュモデルが示す前記複数の節点それぞれの位置座標に基づいて、前記複数の節点の間の距離を示す距離データを生成し、
   前記第１のメッシュモデルの中の複数通りの節点の組それぞれについて、前記距離データが示す前記節点の組に対応する距離に依存する閾値を算出し、前記閾値と前記第１の行列データが示す前記節点の組に対応する係数値とを比較し、比較結果に応じて前記複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新し、
   前記複数の評価値に基づいて前記複数の節点のうちの一部の節点を選択し、前記第１のメッシュモデルから前記選択した一部の節点を除外した第２のメッシュモデルについての第２の行列データを生成し、
   前記第１の行列データおよび前記第２の行列データを用いて、前記複数の節点に対応付けられた複数の変数の解を算出する、
   シミュレーション方法。
5. 複数の節点を含む第１のメッシュモデルについて、前記複数の節点の間の特性を示す係数値を成分として含む第１の行列データを記憶する記憶部と、
   前記第１のメッシュモデルが示す前記複数の節点それぞれの位置座標に基づいて、前記複数の節点の間の距離を示す距離データを生成し、前記第１のメッシュモデルの中の複数通りの節点の組それぞれについて、前記距離データが示す前記節点の組に対応する距離に依存する閾値を算出し、前記閾値と前記第１の行列データが示す前記節点の組に対応する係数値とを比較し、比較結果に応じて前記複数の節点に対応する複数の評価値のうちの少なくとも一部の評価値を更新し、前記複数の評価値に基づいて前記複数の節点のうちの一部の節点を選択し、前記第１のメッシュモデルから前記選択した一部の節点を除外した第２のメッシュモデルについての第２の行列データを生成し、前記第１の行列データおよび前記第２の行列データを用いて、前記複数の節点に対応付けられた複数の変数の解を算出する処理部と、
   を有する情報処理装置。

Images