(57)【要約】【課題】  境界要素法では難しい界面の変化に対してはより柔軟に取り扱うことができるとともに、有限要素法で必要とされる要素分割の数を大幅に減らすことができ、その結果として計算コストを大幅に低減することができる情報処理装置及び情報処理方法並びにプログラムを提供する。【解決手段】  メモリ６２上に、数値解析の対象となる領域を設定し、当該領域内に同一次元の部分領域を設定する。次に、部分領域を要素分割して離散化された基準点を設定する。次に、領域上の物理量と部分領域上の物理量との関係を定義する偏微分方程式Ｌｕ＝ｑを基本解を用いて離散化し、連立一次方程式に変換し、部分領域内の各基準点上で定義される物理量Ｑiを算出する。算出された部分領域内の物理量を領域に対して積分することにより、領域内における部分領域の補集合に相当する領域の物理量ｕを算出する。【選択図】  図２
2005258812.tif 000002
【特許請求の範囲】
【請求項１】
領域を要素分割して数値解析を行う情報処理装置であって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定手段と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手段と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手段と、
前記領域上の物理量と前記部分領域上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する変換手段と、
前記連立一次方程式に基づいて前記部分領域内の各基準点上で定義される物理量を算出する第１の算出手段と、
算出された前記物理量を前記基準点に対応付けて記憶する記憶手段と、
前記記憶手段に記憶された前記部分領域内の物理量を前記領域に対して積分することにより、前記領域内における前記部分領域の補集合に相当する領域の物理量を算出する第２の算出手段と、
前記第２の算出手段による算出結果を出力する出力手段と
を備えることを特徴とする情報処理装置。
【請求項２】
前記物理量が、２次元又は３次元のベクトル量であることを特徴とする請求項１に記載の情報処理装置。
【請求項３】
前記第３の設定手段が、
前記領域をＭ×Ｎに分割する分割手段と、
前記領域内における（Ｍ＋１）×（Ｎ＋１）点分の格子点のうち、前記部分領域に属する格子点を前記基準点として選択する選択手段と
を備えることを特徴とする請求項１又は２に記載の情報処理装置。
【請求項４】
前記変換手段は、
前記領域上で定義される物理量をｕ、線形偏微分作用素をＬ、前記部分領域上の物理量をｑ、前記部分領域に属する前記格子点上の物理量をＱ、前記物理量ｕの近似解をＵとした場合、Ｌｕ＝ｑで表される前記偏微分方程式を離散化して変換されたΣＫＱ＝Ｕで表される前記連立一次方程式におけるマトリックスＫの要素Ｋijは、ｘ、ｙを前記部分領域内の座標として、前記線形偏微分作用素Ｌと２点デルタ関数δ（ｘ−ｙ）から成る方程式Ｌｕ*＋δ（ｘ−ｙ）＝０の基本解ｕ*（ｘ，ｙ）、前記部分領域上の関数の集合ｆi（ｘ）及びｇj（ｙ）を用いて、
Ｋij＝∫ｕ*（ｘ，ｙ）ｆi（ｘ）ｇj（ｙ）ｄｘｄｙ
で定義されることを特徴とする請求項３に記載の情報処理装置。
【請求項５】
前記マトリックスＫは、大部分の要素が零でない値を有する行列要素からなることを特徴とする請求項４に記載の情報処理装置。
【請求項６】
前記線形偏微分作用素Ｌが、ラプラス演算子Δであることを特徴とする請求項４又は５に記載の情報処理装置。
【請求項７】
前記部分領域内に属する前記格子点に通し番号を付与する番号付与手段をさらに備え、
前記第１の算出手段は、前記通し番号が付与された前記部分領域内に属する格子点上の物理量の近似解を算出する
ことを特徴とする請求項３に記載の情報処理装置。
【請求項８】
前記分割手段が、前記領域をＭ×Ｎの正方格子に分割することを特徴とする請求項３に記載の情報処理装置。
【請求項９】
前記部分領域内に設定された基準点に基づいて、該基準点と異なる階層の基準点をさらに細かく設定する第４の設定手段をさらに備えることを特徴とする請求項１から８までのいずれか１項に記載の情報処理装置。
【請求項１０】
前記第４の設定手段が、階層の異なる前記基準点をレベルセット法を用いて設定することを特徴とする請求項９に記載の情報処理装置。
【請求項１１】
コンピュータに領域を要素分割して数値解析を行わせる情報処理方法であって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定工程と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定工程と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定工程と、
前記領域上の物理量と前記部分領域上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する変換工程と、
前記連立一次方程式に基づいて前記部分領域内の各基準点上で定義される物理量を算出する第１の算出工程と、
算出された前記物理量を前記基準点に対応付けて所定の記憶装置に記憶する記憶制御工程と、
前記記憶装置に記憶された前記部分領域内の物理量を前記領域に対して積分することにより、前記領域内における前記部分領域の補集合に相当する領域の物理量を算出する第２の算出工程と、
前記第２の算出工程による算出結果を出力する出力工程と
を有することを特徴とする情報処理方法。
【請求項１２】
コンピュータに、領域を要素分割して数値解析を行わせるためのプログラムであって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定手順と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手順と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手順と、
前記領域上の物理量と前記部分領域上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する変換手順と、
前記連立一次方程式に基づいて前記部分領域内の各基準点上で定義される物理量を算出する第１の算出手順と、
算出された前記物理量を前記基準点に対応付けて所定の記憶装置に記憶させる記憶制御手順と、
前記記憶装置に記憶された前記部分領域内の物理量を前記領域に対して積分することにより、前記領域内における前記部分領域の補集合に相当する領域の物理量を算出する第２の算出手順と、
前記第２の算出手順による算出結果を出力する出力手順と
を実行させるためのプログラム。
【請求項１３】
請求項１２に記載のプログラムを格納したことを特徴とするコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、微分方程式の数値解法において要素分割を用いて近似解を求める情報処理技術に関する。
【背景技術】
【０００２】
微分方程式を数値的に解く手法の中心技術として、有限要素法（ＦＥＭ）や境界要素法（ＢＥＭ）が良く知られている。
【０００３】
ＦＥＭは、考察領域全体を有限要素に分割し、重み付け残差法等を用いて、微分方程式を有限要素式に変換する。次に、当該有限要素式の中の未知数を形状関数と離散化された関数で近似し、要素節点での関数値に関する連立一次方程式を得る。そして、この連立一次方程式を解くことによって、領域内節点での値を求めることができる（例えば、非特許文献１参照）。
【０００４】
一方、ＢＥＭは、支配微分方程式をこれと等価な境界積分方程式に変換した後、考察領域の境界表面を要素分割することで離散化し、ＦＥＭと同様に連立一次方程式に帰着させて解く（例えば、非特許文献２参照）。
【０００５】
上記ＦＥＭは、与えられた汎関数に対する「変分＝０」の形の原理である変分原理による汎関数の最小化条件から定式化される。汎関数は、各要素内で定義される要素汎関数の和として与えられる。このため、ＦＥＭは、領域全体の場の平均を扱うことで、非常にロバストで強力な数値解法となっている。また、ＦＥＭは、区分的近似関数を採用するため、近接相互作用のみを考慮すればよい。その結果、連立一次方程式の係数行列はきわめて「疎な状態」となり、その性質を有効に用いた数値解析法が計算機上で実現されている。
【非特許文献１】O.C.ツィエンキーヴィッツ著、吉識雅夫・山田嘉昭  監訳  “マトリックス有限要素法”、培風館
【非特許文献２】C.A.ブレビア著、神谷紀生・田中正隆・田中喜久昭  共訳、“境界要素法入門”、培風館
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、ＦＥＭは、領域全体を有限要素に分割する必要があるため、データ数が膨大となり、入力データを作成するための計算コストがかかるという問題がある。また、取り扱う方程式の元数も大きなものとなるため、特に３次元問題や非線形問題では、より少ない元数の方程式を扱うことが計算機資源の面からみて望ましいと考えられている。
【０００７】
一方、ＢＥＭは、境界表面だけを要素分割し、各節点にデータを入力すればよいので、計算の前処理が容易である。また、ＦＥＭと比較して、次元を一つ下げて取り扱うことができるため、方程式の元数を少なくすることができ、計算コストを下げることができる可能性がある。
【０００８】
しかし、ＢＥＭは、境界という太さ（幅）のない（測度零の）領域での積分を扱うことになるので、特異点の問題等が発生し、数値計算上ロバストではない場合がある。また、境界積分方程式に変換するために遠距離相互作用を取り入れねばならず、離散化して得られる連立一次方程式の係数行列は「密な状態」となる。その結果、ＢＥＭは、元数は少なくて済む一方で、一般にフルマトリックスを取り扱う必要があり、数値計算上のネックと考えられている。
【０００９】
さらに、ＢＥＭには次のような問題点が存在する。例えば、図１１は、界面の変形により要素の大きさが変化する例を説明するための図である。すなわち、図１１に示すように、粒子が成長して界面の形状が変化するような問題である。この場合、境界表面に節点を置くため、粒子が成長すると最初に設定した節点間の長さが変化し、要素の大きさが変化してしまう。その結果、計算の実行に従って、解析精度が変化してしまうことになる。この問題を避けるためには、要素の大きさが一定になるように、節点の数と配置を計算の実行に合わせて調整する必要があり、そのためには新たな節点へのデータの受け渡しや要素数の変化に対応した行列要素の再計算等をする必要が生じてしまう。
【００１０】
また、図１２は、界面の融合により節点の数が変化する例を説明するための図である。すなわち、図１２に示すように、２つの領域がぶつかって一つの領域に合体するような現象を解析しようとすると、融合した近傍で節点を消滅させたり、新たに発生させたりしなければならない。このように、領域の大きさが変化したり、複数の領域が合体したりするような界面の変化を伴う現象を連続して計算しようとする場合、ＢＥＭでは特有の計算手法が必要となる。
【００１１】
そこで、上述したＦＥＭやＢＥＭのそれぞれの問題を回避するために、ＦＥＭとＢＥＭをハイブリッドに利用するアルゴリズムがこれまでにも提案・研究されている。しかしながら、それらの方法には、特異性の問題等の様々な問題を含んでいる。例えば、Qiya Hu and Dehao Yu, "Solving singularity problems in unbounded domains by coupling of natural BEM and composite grid FEM"  Applied Numerical Mathematics 37(2001)127等がある。
【００１２】
本発明は、このような事情を考慮してなされたものであり、境界要素法では難しい界面の変化に対してはより柔軟に取り扱うことができるとともに、有限要素法で必要とされる要素分割の数を大幅に減らすことができ、その結果として計算コストを大幅に低減することができる情報処理装置及び情報処理方法並びにプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１３】
上記課題を解決するために、本発明は、領域を要素分割して数値解析を行う情報処理装置であって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定手段と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手段と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手段と、
前記領域上の物理量と前記部分領域上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する変換手段と、
前記連立一次方程式に基づいて前記部分領域内の各基準点上で定義される物理量を算出する第１の算出手段と、
算出された前記物理量を前記基準点に対応付けて記憶する記憶手段と、
前記記憶手段に記憶された前記部分領域内の物理量を前記領域に対して積分することにより、前記領域内における前記部分領域の補集合に相当する領域の物理量を算出する第２の算出手段と、
前記第２の算出手段による算出結果を出力する出力手段と
を備えることを特徴とする。
【００１４】
また、本発明は、コンピュータに領域を要素分割して数値解析を行わせる情報処理方法であって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定工程と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定工程と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定工程と、
前記領域上の物理量と前記部分領域上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する変換工程と、
前記連立一次方程式に基づいて前記部分領域内の各基準点上で定義される物理量を算出する第１の算出工程と、
算出された前記物理量を前記基準点に対応付けて所定の記憶装置に記憶する記憶制御工程と、
前記記憶装置に記憶された前記部分領域内の物理量を前記領域に対して積分することにより、前記領域内における前記部分領域の補集合に相当する領域の物理量を算出する第２の算出工程と、
前記第２の算出工程による算出結果を出力する出力工程と
を有することを特徴とする。
【００１５】
さらに、本発明は、コンピュータに、領域を要素分割して数値解析を行わせるためのプログラムであって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定手順と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手順と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手順と、
前記領域上の物理量と前記部分領域上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する変換手順と、
前記連立一次方程式に基づいて前記部分領域内の各基準点上で定義される物理量を算出する第１の算出手順と、
算出された前記物理量を前記基準点に対応付けて所定の記憶装置に記憶させる記憶制御手順と、
前記記憶装置に記憶された前記部分領域内の物理量を前記領域に対して積分することにより、前記領域内における前記部分領域の補集合に相当する領域の物理量を算出する第２の算出手順と、
前記第２の算出手順による算出結果を出力する出力手順と
を実行させるためのプログラムであることを特徴とする。
【発明の効果】
【００１６】
本発明によれば、境界要素法では難しい界面の変化に対してはより柔軟に取り扱うことができるとともに、有限要素法で必要とされる要素分割の数を大幅に減らすことができ、その結果として計算コストを大幅に低減することができる。また、境界要素法では難しかった界面の変化を伴う現象の解析に対して、幅を有する界面として部分領域を設定するとともに、階層的な解像度格子を用いることにより、形状をより細かく表現して、その領域の物理量を連続的に変化させることで、数値計算の取り扱いも容易になるという効果が得られる。
【発明を実施するための最良の形態】
【００１７】
＜第１の実施形態＞
［情報処理装置の構成］
以下、図面を参照して、本発明の実施形態に係る情報処理装置及び当該装置を用いた情報処理方法についてより具体的に説明する。図６は、本発明の第１の実施形態に係る微分方程式の数値解法において要素分割を用いて近似解を求める情報処理装置の構成を示すブロック図である。図６に示すように、本実施形態に係る情報処理装置は、各種演算を実行するＣＰＵ６１やメモリ６２等から構成される演算処理部１、ハードディスク６３やフレキシブルディスク６４からなる記憶部２、キーボード６５やマウス６６等からなる入力部３、そしてディスプレイ６７等からなる出力表示部４からで構成されている。
【００１８】
［情報処理装置の基本的動作］
図２は、本発明の第１の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。本実施形態に係る情報処理装置では、メモリ６２上に図２に示したフローチャートに沿って設計、コード化されたプログラムをロードし、必要な計算領域を確保し、適当な方法で入力された領域形状や格子点座標及び物理量に基づいて所定の演算処理を行う。そして、その結果得られた領域上での物理量の値をハードディスク６３等に書き込んで記憶させたり、ディスプレイ６７上に表示させる。
【００１９】
まず、図１に示すように物理量を演算する対象となる領域Ａを設定する（ステップＳ１）。図１は、本発明の第１の実施形態を説明するための領域Ａとその部分領域ａとを説明するための概念図である。尚、本実施形態に係る領域Ａの次元は２次元としているが、本発明の適用は２次元に限られるものではなく、３次元であってもよい。
【００２０】
次に、領域Ａの部分領域ａを設定する（ステップＳ２）。ここで、部分領域ａは、領域Ａと同一次元の領域である。前述したように、ＦＥＭを用いて領域Ａの物理量を計算する場合は、通常、領域Ａ内の全格子点を未知数として計算しなければならない。また、ＢＥＭを用いて領域Ａの物理量を計算する場合は、部分領域ａに相当するものは次元が一つ少ないものを用いる。しかし、本実施形態では、領域Ａが２次元の場合は部分領域ａは２次元（面）であり、領域Ａが３次元（（立体）の場合は部分領域ａも立体であるとする。また、部分領域ａの形状は任意である。
【００２１】
そして、ｕを領域Ａ上で定義された物理量を表す関数、Ｌを線形偏微分作用素としたとき、領域Ａ上で定義される支配偏微分方程式が、部分領域ａ上で定義されたｑに対して、次式を満たすｑを求めることを考える。
【００２２】
【数１】
2005258812.tif 000003
【００２３】
そして、上記偏微分方程式の解を近似的に求めるため、部分領域ａ上で離散化された格子点を基準点として設定する。尚、各格子点には番号（ｉ＝１，２，…，Ｎ）を付ける（ステップＳ３）。ここで、Ｎは、部分領域ａ上に設定された格子点（例えば、図１において黒丸で示されている格子点）の総数である。尚、基準点となる格子点の離散化を正方格子ではなく、三角格子や任意の非正方規則格子を用いて行っても同様に処理することができる。
次に、物理量ｕを近似表現した各格子点上の値Ｕi（ｉ＝１，２，…，Ｎ）を設定する（ステップＳ４）。尚、格子点上だけに限られるものではなく、格子辺上で定義してもよいし、また格子胞上等で定義してもよい。ここで、次式で示す線形偏微分作用素Ｌと、２点デルタ関数δ（ｘ−ｙ）からなる方程式は、基本解として２点関数ｕ*（ｘ，ｙ）を持つことが知られている。
【００２４】
【数２】
2005258812.tif 000004
【００２５】
ここで、ｘ及びｙは、部分領域ａ内の座標を表す。
また、線形偏微分作用素Ｌとして、ラプラス型やヘルムホルツ型がよく知られているが、本実施形態においてはこれらに限られるものではない。ラプラス型及びヘルムホルツ型には、それぞれに対応して基本解ｕ*（ｘ，ｙ）が求められている。尚、基本解については、前述した非特許文献１に記されている。
【００２６】
そして、上記基本解を用いて、式（１）に示される偏微分方程式を離散化して、次に示すような連立一次方程式に変換する（ステップＳ５）。
【００２７】
【数３】
2005258812.tif 000005
【００２８】
ここで、Ｑj（ｊ＝１，２，…，Ｎ）は、部分領域ａ上の各格子点での値であり、｛Ｑj｝データでｑを近似することができる。また、Ｋij（ｉ，ｊ＝１，２，…，Ｎ）について、基本解ｕ*（ｘ，ｙ）を積分核として次のように定義する。
【００２９】
【数４】
2005258812.tif 000006
【００３０】
ここで、ｆi（ｘ）及びｇj（ｙ）は、それぞれ部分領域ａ上に定義された非零の領域が格子点ｉ又はｊの周りに局在した関数であり、積分は部分領域ａ上で計算される。また、上記積分は解析的又は数値的に行ってもよく、数値積分を実行する際は、ｆi（ｘ）及びｇj（ｙ）を適当な内挿関数とする。尚、このマトリックスＫijは、その大部分の要素（例えば、９０パーセント以上の要素）が非零の値を持つ行列要素からなることを特徴とする。
そして、式（３）で与えられた連立方程式を解くことにより、部分領域ａ上の各格子点上で定義された物理量Ｑiを求めることができる（ステップＳ６）。尚、上述した手順は、部分領域ａ上のＮ個の格子上の値Ｕiから、当該格子上のＱiを求める手続きを示したものであるが、逆にＱiを与えてＵiを求めても良いことは明らかである。
【００３１】
このようにして求められた部分領域ａ上の格子上の値Ｑiから、領域Ａにおける部分領域ａの補集合上であって、座標ｘPを有する任意の点Ｐ上での物理量ｕ（ｘP）は、次式で近似される。
【００３２】
【数５】
2005258812.tif 000007
【００３３】
尚、和は、部分領域ａ上のすべてのＮ個の格子点にわたってとるものとする。
【００３４】
【数６】
2005258812.tif 000008
【００３５】
また、ｇj（ｙ）は前述した関数と同様であり、積分は部分領域ａ上で行うものとする。
以上の手順により、領域Ａの部分領域であるａの離散化された格子点に設定された値から、部分領域ａの領域Ａでの補集合上の任意の点における物理量ｕを求めることができる（ステップＳ７）。
【００３６】
上述したように、通常のＢＥＭによる計算では、部分領域ａとして領域Ａより１つ次元の少ない「表面」を設定し、その上に離散化された節点に物理量を与えるのに対し、本実施形態では、領域Ａと同一次元を有したまま物理量を求めている。従って、ＢＥＭ（境界要素法）では難しかった界面の変化をより柔軟に取り扱うことができる。すなわち、幅を有する界面として部分領域を設定して、その領域の物理量を連続的に変化させることで、数値計算の取り扱いも容易になる。また、領域Ａ内の全格子点について未知数を求めるのではなく、部分領域ａ内の格子点について未知数を設定するので、ＦＥＭ（有限要素法）で必要とされる要素分割の数を大幅に減らすことができ、これによって有限要素法における連立方程式の数を大幅に減らして、計算コストを大幅に低減することができる。
【００３７】
＜第２の実施形態＞
次に、上述した第１の実施形態に係る情報処理装置を用いて、例えば、電磁場のような物理場の物理量を求めるプログラムの流れについて説明する。図７は、本発明の第２の実施形態において計算対象となる領域Ｂ及びその部分領域ｂを説明するための図である。本実施形態では、図７に示すように、領域Ｂを２次元で正方な領域とし、部分領域ｂを同心円で囲まれた領域とした。尚、領域Ｂ及び部分領域ｂの形状は、これだけに限定されることはなく、任意の次元の領域を設定することが可能である。尚、本実施形態における情報処理装置は、前述した第１の実施形態における情報処理装置と同一であり、その処理手順を示すフローチャートは、前述した第１の実施形態で用いられた図２のフローチャートと同一である。すなわち、ステップＳ１、Ｓ２において、図７に示す領域Ｂ及びその部分領域Ｂが設定される。
【００３８】
まず、正方領域ＢをＭ×Ｍの等間隔に分割し、メモリ６２上に（Ｍ＋１）×（Ｍ＋１）点分の格納スペースを作り、そこに各格子点の座標を収納する。また、ｕやｑの離散版であるＵやＱを（Ｍ＋１）×（Ｍ＋１）ベクトルとしてメモリ６２上にとる。次に、各格子点についてその座標が部分領域ｂ上に存在するかどうかの判定を行い、存在する場合のみ選択して、新たなメモリ領域に登録する。その結果、部分領域ｂ上に存在するＮ個の格子点座標ｘiからなる、メモリ空間で連続したベクトルＸを得ることができる。上記処理は、前述のステップＳ３に相当する。尚、正方領域ではない領域を分割する場合は、当該領域をＭ×Ｎに分割し、ＵやＱを（Ｍ＋１）×（Ｎ＋１）ベクトルとしてメモリ上にとればよい。
【００３９】
次に、メモリ空間にＮ個のデータを格納するリストベクトルＬ（ｉ）を用意し、リストベクトルＬ（ｉ）が（Ｍ＋１）×（Ｍ＋１）格子内における部分領域ｂのＮ個のデータを指すようにする。以下、簡単のために、Ｕ（Ｌ（ｉ））をＵi、Ｑ（Ｌ（ｉ））をＱiと記す。但し、ｉ＝１，２，・・・，Ｎである。また同様に、行列要素Ｋijを格納するためのマトリックスＫをメモリ６２上に用意する。
【００４０】
図８は、第２の実施形態において部分領域ｂの格子点上に設定した物理量Ｕiを示す鳥瞰図である。図８では、図を見やすくするために、部分領域ｂ外に対しては一定の値が表示されているが、実際に値は設定されていない。すなわち、この物理量の設定は前述のステップＳ４に相当する。
【００４１】
次に、線形微分作用素Ｌとして２次元ラプラス演算子を選択すると、式（２）から基本解として、次式が得られる。
【００４２】
【数７】
2005258812.tif 000009
【００４３】
次に、マトリックＫの成分を求めるため、式（４）における関数ｆi（ｘ）及びｇj（ｙ）を次式のようにおく。
【００４４】
【数８】
2005258812.tif 000010
【００４５】
【数９】
2005258812.tif 000011
【００４６】
尚、本実施形態では、関数ｆi（ｘ）及びｇj（ｙ）をδ関数としているが、これに限られるものではなく、部分領域ｂで定義された非零の関数であれば良い。上式において、ｘi及びｙjは、それぞれ部分領域ｂ上のｉ番目及びｊ番目の格子点座標である。また、式（７）、式（８）及び式（９）より式（４）を用いて、Ｋijを次のように計算する。
【００４７】
【数１０】
2005258812.tif 000012
【００４８】
すなわち、このＫijの計算は、前述のステップＳ５に相当する。上記計算をするために、ベクトルＸのｉ番目とｊ番目を参照して格子点座標を取り出し、計算された行列要素をマトリックスＫの（ｉ，ｊ）番地に格納する。
次に、ベクトルＵ、マトリックスＫに対して、ベクトル及びマトリック演算の汎用計算ルーチンを使うことにより、式（３）の連立方程式を解いて、部分領域ｂ上の格子点での値Ｑiが格納されたベクトルＱを求める。図９は、第２の実施形態で計算された部分領域ｂ上で定義された物理量ｑを示す鳥瞰図である。図９においても、図８と同様に、部分領域ｂ外に対しては一定の値を入れ、図を見やすくしている。すなわち、このベクトルＱの演算は、前述したステップＳ６に相当する。
【００４９】
次に、部分領域ｂ上の各格子点に対して求められたＱiを用いて、式（５）及び式（６）から、部分領域ｂの領域Ｂでの補集合上の座標ｘPを有する任意の点Ｐに対して物理量ｕ（ｘP）を求め、それを最初に設定した領域ｂ上の格子点上での値Ｕiと合わせて表示する。図１０は、第２の実施形態において領域Ｂ上で計算される物理量ｕを表す鳥瞰図である。すなわち、物理量ｕの取得が、前述したステップＳ７の処理に相当する。
【００５０】
以上のように、計算機上のメモリ空間に設定したベクトルやマトリックスを演算することで、部分領域ｂ上の格子点上の値Ｕiから領域Ｂ上の物理量ｕを近似的に求めることができる。
【００５１】
＜第３の実施形態＞
図３は、本発明の第３の実施形態に係る情報処理装置による処理内容を説明するための図である。本実施形態における情報処理装置の構成及びその動作は前述したものと同様である。図３において、領域Ｃの部分領域ｃは２つの分離された領域からなる。
【００５２】
本実施形態においても、前述した第２の実施形態と同様に部分領域ｃ上で離散化された格子点を設定する。そして、設定された格子点数をＮとして、各格子点上に物理量を設定することにより、第２の実施形態と同様の処理手順に従って、領域Ｃ上の任意の点における物理量ｕを求めることができる。
【００５３】
尚、部分領域ｃが３つ以上の分離された領域からなる場合であっても、上記手順が同様に成り立つことは明らかである。
【００５４】
＜第４の実施形態＞
図４は、本発明の第４の実施形態に係る情報処理装置による処理内容を説明するための図である。本実施形態における情報処理装置の構成及びその動作は前述したものと同様である。図４に示すように、本実施形態では、領域Ｄは、その部分領域ｄの境界の一部が一致している。
【００５５】
本実施形態においても、前述した第２の実施形態と同様に、まず、部分領域ｄ上で離散化された格子点を設定する。そして、設定された格子点数をＮとして、各格子点上に物理量を設定することにより、第２の実施形態と同様の処理手順に従って、領域Ｄ上の任意の点における物理量ｕを求めることができる。
【００５６】
尚、部分領域ｄの境界が領域Ｄの境界全部と一致している場合でも、上記手順が同様に成り立つことは明らかである。
【００５７】
＜第５の実施形態＞
図５は、本発明の第５の実施形態に係る情報処理装置による処理内容を説明するための図である。本実施形態における情報処理装置の構成及びその動作は前述したものと同様である。図５に示すように、本実施形態では、部分領域ｅ上に非正方格子状に格子点を設定した場合を示す。この場合も、前述した第２の実施形態と同様に、まず、部分領域ｅ上の適当な位置にＮ個の離散化された格子点を設定する。そして、このときの各格子点上に物理量を設定することにより、第２の実施形態と同様の処理手順に従って、領域Ｅ上の任意の点における物理量ｕを求めることができる。
【００５８】
また、格子点の離散化を正方格子ではなく、三角格子や任意の非正方規則格子を用いて行っても同様に処理することができる。さらに、適当な擬似乱数を用いて格子点位置を設定した場合であっても、上記手順が同様に成り立つことは明らかである。
【００５９】
＜第６の実施形態＞
前述した各実施形態では、部分領域上の格子上に物理量Ｕを与えて、Ｑを解いたが、部分領域上のＱの一部を既知量として、残りのＱを解くようにしても良い。このとき、式（３）で与えられる連立方程式の数と未知量の数が一致するように、Ｕの一部を未知量として取り扱う。
【００６０】
＜第７の実施形態＞
前述した第１の実施形態では、図２に示すように、まず領域Ａを設定し（ステップＳ１）、次に領域Ａと同一次元の部分領域ａを設定した（ステップＳ２）。本実施形態では、部分領域ａの形状をより良く表現して計算する。すなわち、本実施形態では、部分領域ａ内の格子を階層的な多重格子を利用して、格子をさらに細かくする。
【００６１】
このような、所望の精度で必要な領域の解像度を上げる計算方法としては、例えば、Jhon  Strain,  " Fast Tree-Based Redistancing for Level Set Computations", Journal of Computational Physics 152, p.664-686(1999)や、Michael Griebel, Gerhard Zumbusch, "Parallel adaptive subspace correction schemes with applications to elasticity" Computer Methods in applied mechanics and engineering 184 p.303-332 (2000)に提案されている。そこで、本実施形態では、これらの方法を用いて実現する。
【００６２】
図１３は、本発明の第７の実施形態を説明するための領域Ｆとその部分領域ｆとを説明するための概念図である。また、本実施形態における情報処理装置の構成及びその動作は第１の実施形態で用いられた情報処理装置と同様である。また、図１４は、本発明の第７の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。ここで、領域Ｆの設定（ステップＳ２１）は第１の実施形態でのステップＳ１と、部分領域ｆの設定（ステップＳ２２）はステップＳ２の処理と同様である。
【００６３】
第７の実施形態では、部分領域ｆを設定した後、領域Ｆとは異なる階層で、さらにその部分領域ｆ内の格子の細分化を行う（ステップＳ２３）。そして、細分化された部分領域ｆ上で離散化した格子点の設定／番号付け処理（ステップＳ２４）、番号付けされた格子点への物理量Ｕiの割り当て処理（ステップＳ２５）、マトリックスＫの作成処理（ステップＳ２６）、部分領域ｆ上の格子点上で定義される物理量Ｑiの算出処理（ステップＳ２７）、領域Ｆ上の物理量ｕの算出処理（ステップＳ２８）は、第１の実施形態のステップＳ３〜Ｓ７の処理を細分化された格子点に対して同様に行うものである。
【００６４】
このように、本実施形態では、部分領域ｆ内の格子を第１の実施形態の場合よりも細かく設定したので、領域Ｆの元の格子よりも小さな変化を有する形状に対しても対応することが可能である。但し、部分領域ｆ内の解像度が異なる場合は、胞内の体積の考慮が必要である。
【００６５】
＜第８の実施形態＞
次に、上述した第７の実施形態に係る情報処理装置を用いて、例えば、電磁場のような物理場の物理量を求めるプログラムの流れについて説明する。本実施形態で用いる領域と部分領域は、前述した第２の実施形態で用いられた図７に示す領域Ｂと同心円で囲まれた部分領域ｂを用いる。尚、領域Ｂ及び部分領域ｂの形状は、これだけに限定されることはなく、任意の次元の領域を設定することが可能である。尚、本実施形態における情報処理装置は、前述した第１、第７の実施形態における情報処理装置と同一であり、その処理手順を示すフローチャートは、前述した第７の実施形態で用いられた図１４のフローチャートと同一である。すなわち、ステップＳ２１、Ｓ２２において、図７に示す領域Ｂ及びその部分領域Ｂが設定される。
【００６６】
そして、正方領域ＢをＲ×Ｒの等間隔に分割し、さらに部分領域を有する格子を細分化する（ステップＳ２３）。本実施形態では、最上位階層を１とすると、２階層まで細分化している。尚、この階層数は、この限りでなく、任意の階層数を設定することができる。そして、このように多重化した格子を含めて、メモリ６２上にＭ×Ｍ点分の格納スペースを作り、そこに各格子点の座標を収納する。また、ｕやｑの離散版であるＵやＱをＭ×Ｍベクトルとしてメモリ６２上にとる。次に、各格子点についてその座標が部分領域ｂ上に存在するかどうかの判定を行い、存在する場合のみ新たなメモリ領域に登録する。その結果、部分領域ｂ上に存在するＮ個の格子点座標ｘiからなる、メモリ空間で連続したベクトルＸを得ることができる。上記処理は、前述のステップＳ２４に相当する。尚、正方領域ではない領域を分割する場合は、当該領域をＭ×Ｎに分割し、ＵやＱを（Ｍ＋１）×（Ｎ＋１）ベクトルとしてメモリ上にとる。
【００６７】
次に、メモリ空間にＮ個のデータを格納するリストベクトルＬ（ｉ）を用意し、リストベクトルＬ（ｉ）がＭ×Ｍ格子内における部分領域ｂのＮ個のデータを指すようにする。以下、簡単のために、Ｕ（Ｌ（ｉ））をＵi、Ｑ（Ｌ（ｉ））をＱiと記す。但し、ｉ＝１，２，・・・，Ｎである。また同様に、行列要素Ｋijを格納するためのマトリックスＫをメモリ６２上に用意する。
【００６８】
以下、前述した第２の実施形態と同様に、領域Ｂ上の物理量と部分領域ｂ上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化し、連立一次方程式に変換する。そして、連立一次方程式に基づいて部分領域ｂ内の各格子点上で定義される物理量を算出し、部分領域ｂ内の物理量を領域Ｂに対して積分することにより、領域Ｂ内における部分領域ｂの補集合に相当する領域の物理量を算出することができる。
【００６９】
以上のように、本実施形態では、計算機上のメモリ空間に設定したベクトルやマトリックスを演算することで、部分領域ｂ上のより細かく設定された格子点上の値Ｕi から領域Ｂ上の物理量ｕを近似的に求めることができる。
【００７０】
また、前述した第３、第４の実施形態と同様に、部分領域が２つ以上の領域から成る場合であっても、領域と部分領域の一部の境界が一致している場合であっても、本発明を同様に適用することが可能である。すなわち、図１５は、領域Ｇ内に２つ以上の部分領域ｇが存在する場合の細分化された格子の状態を説明するための図である。また、図１６は、領域Ｈと部分領域ｈの一部に境界線が一致する場合の細分化された格子の状態を説明するための図である。
【００７１】
＜第９の実施形態＞
図１７及び図１８は、本発明の第９の実施形態に係る情報処理装置による処理内容を説明するための図である。図１７に示すように、本発明の多くの問題においては、部分領域は何らかの界面Γ（領域Ｒが３次元の場合の界面Γは２次元曲面、領域Ｒが２次元の場合の界面Γは曲線）を所望の太さ（幅）を持たせるようにする膨張工程を含んだものとなる。
【００７２】
そこで、本実施形態では、「Jhon Strain, " Fast Tree-Based Redistancing for Level Set Computations", Journal of Computational Physics 152, p.664-686(1999)」のレベルセット法を用いて、図１７に示すように階層格子を生成するものとする。具体的には、界面Γからの距離関数をレベルセット法によって作成し、この関数に基づいて格子を分割する。図１７に示す例では、最上位階層を１とすると、４階層まで細分化している。そして、界面生成後は、図１８に示すように、階層格子の解像度を考慮し、領域を選択する。尚、図１８に示す例では、第３階層以上（３階層と４階層を含む領域）を持つ部分領域をｒとしている。
【００７３】
また、階層格子を作る方法は、本実施形態で示した手法に限られるわけではな。本実施形態の本質は、階層格子により部分領域ｒを作成し、これを利用して本目的の演算を実現することにある。
【００７４】
さらに、本発明では、格子点の離散化を正方格子で行なったが、三角格子や任意の構造格子を用いていってもよい。すなわち、領域内に存在する任意構造の格子に対して、図１４のフローチャートにおけるステップＳ２３において、部分領域を設定した後にその格子の細分化を行ったり、界面Γがある場合は、界面Γからの距離関数を用いて階層格子を作成するようにしてもよい。
【００７５】
＜その他の実施形態＞
以上、実施形態例を詳述したが、本発明は、例えば、システム、装置、方法、プログラム若しくは記憶媒体等としての実施態様をとることが可能であり、具体的には、複数の機器から構成されるシステムに適用しても良いし、また、一つの機器からなる装置に適用しても良い。
【００７６】
尚、本発明は、前述した実施形態の機能を実現するソフトウェアのプログラム（実施形態では図に示すフローチャートに対応したプログラム）を、システムあるいは装置に直接あるいは遠隔から供給し、そのシステムあるいは装置のコンピュータが該供給されたプログラムコードを読み出して実行することによっても達成される場合を含む。
【００７７】
従って、本発明の機能処理をコンピュータで実現するために、該コンピュータにインストールされるプログラムコード自体も本発明を実現するものである。つまり、本発明は、本発明の機能処理を実現するためのコンピュータプログラム自体も含まれる。
【００７８】
その場合、プログラムの機能を有していれば、オブジェクトコード、インタプリタにより実行されるプログラム、ＯＳに供給するスクリプトデータ等の形態であっても良い。
【００７９】
プログラムを供給するための記録媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光ディスク、光磁気ディスク、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、磁気テープ、不揮発性のメモリカード、ＲＯＭ、ＤＶＤ（ＤＶＤ−ＲＯＭ，ＤＶＤ−Ｒ）などがある。
【００８０】
その他、プログラムの供給方法としては、クライアントコンピュータのブラウザを用いてインターネットのホームページに接続し、該ホームページから本発明のコンピュータプログラムそのもの、もしくは圧縮され自動インストール機能を含むファイルをハードディスク等の記録媒体にダウンロードすることによっても供給できる。また、本発明のプログラムを構成するプログラムコードを複数のファイルに分割し、それぞれのファイルを異なるホームページからダウンロードすることによっても実現可能である。つまり、本発明の機能処理をコンピュータで実現するためのプログラムファイルを複数のユーザに対してダウンロードさせるＷＷＷサーバも、本発明に含まれるものである。
【００８１】
また、本発明のプログラムを暗号化してＣＤ−ＲＯＭ等の記憶媒体に格納してユーザに配布し、所定の条件をクリアしたユーザに対し、インターネットを介してホームページから暗号化を解く鍵情報をダウンロードさせ、その鍵情報を使用することにより暗号化されたプログラムを実行してコンピュータにインストールさせて実現することも可能である。
【００８２】
また、コンピュータが、読み出したプログラムを実行することによって、前述した実施形態の機能が実現される他、そのプログラムの指示に基づき、コンピュータ上で稼動しているＯＳなどが、実際の処理の一部または全部を行い、その処理によっても前述した実施形態の機能が実現され得る。
【００８３】
さらに、記録媒体から読み出されたプログラムが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によっても前述した実施形態の機能が実現される。
【図面の簡単な説明】
【００８４】
【図１】本発明の第１の実施形態を説明するための領域Ａとその部分領域ａとを説明するための概念図である。
【図２】本発明の第１の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。
【図３】本発明の第３の実施形態に係る情報処理装置による処理内容を説明するための図である。
【図４】本発明の第４の実施形態に係る情報処理装置による処理内容を説明するための図である。
【図５】本発明の第５の実施形態に係る情報処理装置による処理内容を説明するための図である。
【図６】本発明の第１の実施形態に係る微分方程式の数値解法において要素分割を用いて近似解を求める情報処理装置の構成を示すブロック図である。
【図７】本発明の第２の実施形態において計算対象となる領域Ｂ及びその部分領域ｂを説明するための図である。
【図８】第２の実施形態において部分領域ｂの格子点上に設定した物理量Ｕiを示す鳥瞰図である。
【図９】第２の実施形態で計算された部分領域ｂ上で定義された物理量ｑを示す鳥瞰図である
【図１０】第２の実施形態において領域Ｂ上で計算される物理量ｕを表す鳥瞰図である。
【図１１】界面の変形により要素の大きさが変化する例を説明するための図である。
【図１２】界面の融合により節点の数が変化する例を説明するための図である。
【図１３】本発明の第７の実施形態を説明するための領域Ｆとその部分領域ｆとを説明するための概念図である。
【図１４】本発明の第７の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。
【図１５】領域Ｇ内に２つ以上の部分領域ｇが存在する場合の細分化された格子の状態を説明するための図である。
【図１６】領域Ｈと部分領域ｈの一部に境界線が一致する場合の細分化された格子の状態を説明するための図である。
【図１７】本発明の第５の実施形態に係る情報処理装置による処理内容を説明するための図である。
【図１８】本発明の第５の実施形態に係る情報処理装置による処理内容を説明するための図である。
【符号の説明】
【００８５】
１  演算処理部
２  記憶部
３  入力部
４  出力表示部
６１  ＣＰＵ
６２  メモリ
６３  ハードディスク
６４  フレキシブルディスク
６５  キーボード
６６  マウス
６７  ディスプレイ
【図１】
2005258812.tif 000013
【図２】
2005258812.tif 000014
【図３】
2005258812.tif 000015
【図４】
2005258812.tif 000016
【図５】
2005258812.tif 000017
【図６】
2005258812.tif 000018
【図７】
2005258812.tif 000019
【図８】
2005258812.tif 000020
【図９】
2005258812.tif 000021
【図１０】
2005258812.tif 000022
【図１１】
2005258812.tif 000023
【図１２】
2005258812.tif 000024
【図１３】
2005258812.tif 000025
【図１４】
2005258812.tif 000026
【図１５】
2005258812.tif 000027
【図１６】
2005258812.tif 000028
【図１７】
2005258812.tif 000029
【図１８】
2005258812.tif 000030