先ず最初に、FreeFEM++は有限要素法解析のためのソフトウェアらしいです。Octaveとかそういうのの有限要素法バージョンだと思えば良いんだろうか。
そして俺は有限要素法のことなんざこれっぽっちも分かっちゃいません。ごくごく初歩的なこととして、変分原理を使うことと、重み関数というものが登場することと、微分方程式を特殊な形式に書き直して解析することしか知りません。それで有限要素法のためのコーディングをすると回がポチっとなと出て来ることしか分からない。因みに有限要素法は偏微分方程式の境界値問題を解くための方法らしい。
そんな俺でもFreeFEM++があれば有限要素法使えちゃうんですよ。技術の進歩って凄いですね!!開発はフランスのピエール・マリー・キュリー大学*1のグループが中心となってやってて、日本人の方も開発に携わってるらしいです。
FreeFEM++を手に入れる
から自分の環境にあったのを引っ張ってきます。ソースからコンパイルしたいとかマニアックな趣味の方はそうすれば良いかと。普通にバイナリパッケージを引っ張ってきて、インストールマニュアルどおりにやれば良いと思います。
僕の場合はOS Xのppcだったので、そのバイナリを持ってきて、/Applicationに入れて、何だかパスを通すようなことをやった気がする。でもパスを通す云々はCUIでやりたいときだけ必要っぽいので多分不要。
試しに動かす
その後、FreeFEM++を起動します。
すると、FreeFEM++は[ファイル]と、[ヘルプ]だけのシンプルな構成であることが分かります。
そして[file]から、FreeFEM++のバイナリだとかサンプルだとかを置いてあるディレクトリに行きます。そして任意の.edpファイルを開くと実際に例がでてきます。それでちゃんとグラフが表示されればおk。
自分で問題を解いてみる
そしたら次は実際に自分で問題を解いてみます。
有限要素法解析に持ち込むにはそれなりに式の変形をしなければならないらしいです。取り敢えずFreeFEM++はbilinearな形でないと読み込めないらしい。ので、支配方程式を変形します。
今、方程式をLaplace方程式として、△u=0で、境界上である値をとるという風に考えます。
その次に、変分を取るときだかに必要な重み関数というものを導入します。これをvとします。vは境界上で0に収束します。そして、支配方程式にvを掛けます。それで部分積分すると、
になります。そして、境界Cの上でvは0なので、下の式の右辺の第一項は消えるらしいです。ので、結局方程式は、
になります。
また、境界条件は適当に与えます。まあ楕円の上半分で値を持って、下半分では0とかそんな感じかね。そして今、境界は楕円で置くことにします。ということで、以下FreeFEM++での問題の解き方です。
最初に色んな言語でやるように諸々の定数を決めときます。今、楕円のアスペクト比と値を与える領域のパラメータなんぞを与えます。ということで、
real a = 1.0, b = 2.0; real Theta = 0.5*pi;
ここで、文末はセミコロンで区切るのはやっぱり色んな言語と同じです。realは実数の定数であることを示すみたいです。
その次に境界と領域を設定します。楕円なんですが、FreeFEM++では媒介変数を用いてそれらを定義することが出来ます。そして楕円の上半分をG1、下半分をG2とします。あとpiには最初から円周率が与えられてるみたいです。
border G1 (t = 0, Theta) { x = a*cos(t); y = b*sin(t); }
border G2 (t = Theta, 2*pi) { x = a*cos(t); y = b*sin(t); }です。borderで領域を設定するようです。ここで、t=0からThetaまでがG1で、Thetaから2πまでがG2です。二つに分けるのには境界条件を与える都合上のあれがあるので、楕円そのものを与えたい場合は、
border G (t = 0, 2*pi) { x = a*cos(t); y = b*sin(t); }とかやればおkです。
その次にメッシュを切ります。この辺は良く分かりませんが、サンプルコードを真似っこしてたら出来たんで良しとしましょう。
mesh Th = buildmesh (G1(50) + G2(50)); fespace Vh (Th, p1); Vh u, v;
一番上はどうやらメッシュをいくつきるかというのが決められるらしいのが分かります。G1とG2でそれぞれ50点メッシュを切ってるらしいです。その次の行で、境界Thの中の領域をどうやってメッシュで切るかを決めるらしいです。VhだとかThだとかは何だか良く分からないけど、Thは領域らしいです。P1とかいうのは、P0からP9だか位まであって色んな形のメッシュを切ってくれるらしいです。取り敢えずP1で良しとします。そしてその下の行で未知数がuとvであることを示します。ここで、uが方程式の解で、vが重み関数です。
その次に方程式を書き込みます。境界条件もこの方程式を書き込む所で一緒にぶち込みます。この辺は良くあるCFDのあれに似てますね。
solve Laplace (u, v) = int2d (Th) (dx(u)*dx(v) + dy(u)*dy(v)) + on (G1, u = x);
この辺は方程式の変形とか、境界条件がどうだから非同次になるとかの考察をあんましやってないので微妙ですが、形だけ。
まず、solveというので方程式を解きます。で、関数に名前をつけて、後々使いたいときは、Laplace(u,v)とか名前を付けるみたいです。別になくても構わないみたいです。
そして、int2dというのが、微分演算子が一つの項の中に2つあるときにつけるもののようです。そして、微分はdx(hoge)でhogeをxで微分という風にします。ので、一回微分の項の場合は、int1d(Th)(dx(hoge))とするみたいです。非同次項はどこにでもいいから括弧ん中にぶち込んどけば良いみたいです。
その次に境界条件です。
境界条件はonによって指定されます。
そして、G1の上でuがxと同じ値をとるということを上のでは示しています。
最後に図をプロットするのに、
plot (u);
でやります。
これらをまとめたもの、
real theta = 0.5*pi;
real a = 1.0, b = 2.0;
func z = x;
border G1 (t = 0, theta) { x = a*cos(t); y = b*sin(t); }
border G2 (t = theta, 2*pi) { x = a*cos(t); y = b*sin(t); };
mesh Th = buildmesh (G1(50) + G2(50));
fespace Vh (Th, P1);
Vh u, v;
solve Laplace (u,v) = int2d (Th) (dx(u)*dx(v) + dy(u)*dy(v)) + int1d(Th)(v)
+on (G1, u = z) + on (G2, u = 2*z);
plot (u);を適当に保存して、FreeFEM++から開くと実際にLaplace方程式を解いたものがでてきます。epsファイルで保存したいときは、
plot (u, wait =1, ps = "Laplace.eps");
とかするといいみたいです。
矩形の境界
FreeFEM++では矩形の境界も簡単に設定できます。
領域が閉じてれば良いみたいです。なので、2×2の境界上で、境界条件が右上で何らかの値をとって、左下で0の場合というのは、
border a0 (t = 0, 2) { x = t; y = 0; }
border a1 (t = 0, 2) { x = 2; y = t; }
border a2 (t = 0, 2) { x = 2-t; y = 2; }
border a3 (t = 0, 2) { x = 0; y = 2-t; }として、支配方程式の部分を、
solve Laplace (u,v) = int2d (Th) (dx(u)*dx(v) + dy(u)*dy(v)) + int1d(Th)(v)
+on (a1, a2, u = z);とすれば良いみたいです。
*1:Universit? Pierre et Marie Curieらしい。
