FreeFEM++ の segment の使い方

有限要素法のフリーソフトであるFreeFEM++でsegmentというコマンドを使おうとして苦労したのでメモを作っておく．
→ 最終的な結論：buildmeshLというのがあるらしいのでそっちを使う方が良さそう．

segmentコマンドについてのマニュアル：The command segment

segmentは3次元空間での曲線を定義するときに使う．
使い道の一つとして，線積分をしたいときの経路を決める場合に使える．
その簡単な例をFreeFEM++の開発者のHechtさんが書いてくれている．

一番簡単な使い方

一番簡単な使い方は

segment1.edp

load "msh3";
meshL Th = segment(10,[x,0,0]);
plot(Th);

というファイルを作成して
$ FreeFEM++ segment1.edp
を実行すると以下のようなウィンドウで結果が表示される．

解説すると，

1. segmentを使うためにload "msh3"で3次元メッシュモジュールを読み込む．
2. meshL Th = segment(10,[x,0,0]);で，$x=0$から$x=1$まで動く直線を10分割したメッシュを作る．
3. plot(Th);で結果をプロット．

ちなみに，デフォルトの設定が[x,0,0]らしく，

meshL Th = segment(10);

でも同じことが実行される．

注意点

• meshL Th = segment(10,[x,0,0]);においてxはxでないといけない模様．つまり，他の文字yとかzとかに置き換えるとエラーになった．
• マニュアルにはxがどんな値をとるのか書いておらず，色々調べたところ，$x=0$から$x=1$ということが分かった．

関数による曲線の表示1

[x,0,0]の部分をもう少し工夫して複雑な曲線でメッシュを作る．

たとえば，

segment2.edp

load "msh3";
func Fx = cos(2*pi*x);
func Fy = sin(2*pi*x);
meshL Th = segment(10,[Fx,Fy,0]);
plot(Th);

というファイルを作って
$ FreeFEM++ segment2.edp
を実行すると，以下のようなウィンドウが表示される．

分割数が10なのでかなり綺麗な円には見えないが，分割数を増やせば綺麗な円になる．（分割数100の場合）

注意点

• ここでもfunc Fx = cos(2*pi*x);とfunc Fy = sin(2*pi*x);でxをyとかに変えるとエラーになってしまう．エラーを吐く行がmeshL 〜の行なので，関数の定義の問題というよりは，segmentのところが問題なようだ．

関数による曲線の表示2（うまく動かず）

マニュアルにある例

segment3.edp

load "msh3";
real R = 3, r=1;
real h = 0.1; //
int nx = R*2*pi/h;
func torex= (R+r*cos(y*pi*2))*cos(x*pi*2);
func torey= (R+r*cos(y*pi*2))*sin(x*pi*2);
func torez= r*sin(y*pi*2);
meshL Th=segment(nx,[torex,torey,torez],removeduplicate=true) ;
plot(Th);

を実行してみても，ただの円が表示されるだけだった（FreeFEM++ v4.11)．どうもyがずっと0のままみたい．うーん．

ということで2変数での曲線のメッシュ化は諦めた．