conormalというのは、偏微分方程式を a*u-div(C.grad(u))、法線ベクタを n としたとき、(Cgrad(u))*n という意味です。

内部的には、すべての境界条件は

  id(u) * a + dnu(u) * b = c;

関数 a, b, c は preciseが指定されないときは要素内で連続、指定すると不連続になりま す。

Dirichlet-Neumann corner(１辺がディリクレ、他の辺がノイマンの角) ではpreciseが指定されないときは、関数 a, b, c を要素内で連続にし て、角をディリクレ条件で扱い、preciseが指定されたときは、a, b, c を不連続にして、ノイマン条件の辺が正しく処理されるようにします。

preciseが指定されたときは、評価関数に１次でなく２次式を使い、各頂点 の少し内側の３点(Gauss point)を使ってガウス積分を使うようになります。

偏微分方程式は

  pde(u) [+-]op1(u)[*/]exp1 [+-]op2[*/]exp2...=exp3;

  id(u) - u(x,y) の値
  dx(u) - (d/dx)u(x,y) .. x で偏微分
  dy(u) - (d/dy)u(x,y) .. y で偏微分
  laplace(u) - ((d/dx)^2+(d/dy)^2)u(x,y) .. ラプラシアン
  dxx(u) - (d/dx)((d/dx)u(x,y)) .. x で２回偏微分
  dxy(u) - (d/dx)((d/dy)u(x,y)) .. y で偏微分した後 x で偏微分
  dyx(u) - (d/dy)((d/dx)u(x,y)) .. x で偏微分した後 y で偏微分
  dyy(u) - (d/dy)((d/dy)u(x,y)) .. y で２回偏微分

  laplace(u)*v

  laplace(u)*v - (d/dx+d/dy)(v(x,y)*((d/dx+d/dy)f(x,y))

2.10. 作図

  plot(u);
  plot3d(u);

複素変数のplot(f),plot3d(f)はデフォルトで u の実部を描画します ので、虚部を描画する場合は plot(Im(f)) と指定します。

2.11. 結果の保存

下記のようなデータの保存とロード関数が用意されています。

  save('file',f) - file に f を保存
  savemesh('file') - file にメッシュデータを保存
  load('file',f) - file から f をロード
  loadmesh('file') - file からメッシュデータをロード

save() で f が関数の場合は、次の形式のデータが書き込まれます。

  ns
  f[0]
  f[1]
  ....
  f[ns-1]

savemesh() では次の形式のデータが作られます。

  ns nt
  q[0].x q[0].y ib[0]
  ...
  q[ns-1].x q[ns-].y ib[ns-1]
  me[0][0] me[0][1] me[0][2] ngt[0]
  ..
  me[nt-1][0] me[nt-1][1] me[nt-1][2] ngt[nt-1]

メッシュの適応制御などで、連番を付けたファイル名がほしいときは、下記の構文 を使います。

  save('file',f,c)
  savemesh('file',c)
  load('file',f,c)
  loadmesh('file',c)

2.12. 連立方程式

次の方法で連立偏微分方程式が解けますが、コンピュータの負担は極端に増えま すから、単純な繰り返しで解くことができないことを確認してから使ってくださ い。

  solve(u,v) {
   onbdy(..) ...dnu(u)...=... /* u の境界条件 */
   onbdy(..) u=...            /* u の境界条件 */
   pde(u) ...                 /* u の偏微分方程式 */
   onbdy(..)<v=... または ...dnu(v)...> /* v の境界条件 */
   pde(v) ...                 /* v の偏微分方程式 */
  };

solveの構文はスカラの場合と同じで solve(u,v,1) 等も使えます。方程式 の順序は問いません。onbdy()はいくつでも書けます。onbdy()の 構文もスカラの時と同じですが、dnu()は一つしか書けません。

2.13. 弱形式 (Weak formulation)

次の構文で「弱形式」(weak form)による偏微分方程式の記述ができます。

  varsolve(<unknown function list>;<test function list>,<<int>>)
	<<instruction>> : <<expression>>;

<unknown function list> と <test function list> は１つないし ２つの関数をコンマで区切って指定します。<int> は正か負の整数です。 <instruction> は単文か複文、<expression> は実数または複素数を返す 式です。<< >> で囲まれた項は省略可能です。

次の例を見てください。

varsolve(u;w)   /* ディリクレ問題 -laplace(u) =x*y */
begin
  onbdy(1) u = 0;
  f = dx(u)*dx(w) + dy(u)*dy(w)
  g = x*y;
end : intt[f] - intt[g,w];


varsolve(u;w,-1)   /* 同じ (組み立て済みマトリクスを使用) */
begin
  onbdy(1) u = 0;
  f = dx(u)*dx(w) + dy(u)*dy(w)
  g = x*y;
end : intt[f] - intt[g,w];


varsolve(u;w)   /* ノイマン問題 u-laplace(u) = x*y */
begin
  f = dx(u)*dx(w) + dy(u)*dy(w) -x*y;
  g = x;
end : intt[f] + int[u,w] - int(1)[g,w];


varsolve(u,v;w,s) /* Lame 方程式 */
begin
   onbdy(1) u=0;
   onbdy(1) v=0;
   e11 = dx(u);
   e22 = dy(v);
   e12 = 0.5*(dx(v)+dy(u));
   e21 = e12;
   dive = e11 + e22;
   s11w=2*(lambda*dive+2*mu*e11)*dx(w);
   s22s=2*(lambda*dive+2*mu+e22)*dy(s);
   s12s = 2*mu*e12*(dy(w)+dx(s));
   s21w = s12s;
   a = s11w+s22s+s12s+s21w +0.1*s;
end : intt[a];

varsolveはsolveに比べて計算に時間がかかりますが、これは、 dx(u) のような微分がsolveの 1 度に対して、varsolveは 9 回 も必要になるためです。

2.14. その他の関数

以上に加えて、次のような関数があります。

                         _
  scal(f,g) - int(f(x,y)*g(x,y))*dx*dy .. int() は面積分
  wait - 結果を表示してマウスの左釦が押されるのを待つ
  nowait - マウスの左クリックを待たない
  changewait - wait/nowait を反転
  precise - １次でなく２次の評価関数を使う

a:=scal(f,g) は三角分解された領域について、次の積分(f(x,y) と g(x,y) の 共役複素数の積)を実行します。

                _
  a = ∫(f(x,y)*g(x,y))*dx*dy

preciseはFreeFEMに２次の積分公式を使った高精度の計算を行う ように指示します。積分は個々の３角形の少し内側の３(Gause)点で行われ ますから、積分を行う関数は、より多くのメモリ(関数あたり ns が 3*nt、つまり、 9 倍)を消費しますが、これは、偏微分方程式 (PDE) が必要とするメモリに比べれ ばわずかなものです。

preciseをconvectや不連続な非線形システムで使うのは、とても よい考えです。

FreeFEMのplotコマンドは描画後、ユーザがその結果をよく見られる ようにプログラムを一次停止させ、マウスの左釦がクリックされるのを待ちますが、 nowaitを実行すると、一次停止しないようになり、waitが実行される と再び停止するようになります。changewaitはこの２つの状態を反転させま す。また、X11を使ったシステムでは、マウスの右釦をクリックすると、プ ログラムの実行を停止します。 waitは

2.15. 外部関数とのリンク

exec('command')はUnixのshに指定されたcommandを 実行させます。

user(what,f)はFreeFEM内部の C++ のプログラムで、下記の動作を 行います。

  for (int j = 0; j < nquad; j++)
    creal gfemuser(creal what, creal* f, int j);

whatは複数のユーザ定義関数を必要とするとき、どの関数を使うかを if文等で選べるようにするための引数です。

gfemuser()の中では、三角分解(ns,nt,me,q,ng,ngt,area,..)等、全ての グローバル関数が参照可能です。src/fem.cxx にgfemuserの実装例があり ますので、見てください。

2.16. saveall()

下記の構文で、すべてのデータをまとめてファイルに書き込むことができ、ユーザ が自分の解析プログラムを書けるようになっています。

  saveall('fime_name',var_name1,...);

u=p   ディリクレ
c u+dnu(u)=g ノイマン
b u-dx(nuxx dx(u))-dx(nuxy dy(u))-dy(nuyx dx))-dy(nuyy dy(u))
 + a1 dx(u) + a2 dy(u) =f

実際の C++ のルーチンは次のとおりです。

int saveparam(fcts *param, triangulation* t, char *filename, int N)
{
  int k, ns = t->np;
  ofstream file(filename);
  file<<ns<<"    "<<N<<endl;
  for(k=0; k<ns; k++)
  {
                  file << (param)->f[k]<<" " ; file<<"             ";
                  file << (param)->g[k]<<" " ; file<<"             ";
                  file << (param)->p[k]<<" " ; file<<"             ";
                  file << (param)->b[k]<<" " ; file<<"             ";
                  file << (param)->c[k]<<" " ; file<<"             ";
                  file << (param)->a1[k]<<" " ; file<<"            ";
                  file << (param)->a2[k]<<" " ; file<<"            ";
                  file << (param)->nuxx[k]<<" " ; file<<"          ";
                  file << (param)->nuxy[k]<<" " ; file<<"          ";
                  file << (param)->nuyx[k]<<" " ; file<<"          ";
                  file << (param)->nuyy[k]<<" " ; file<<"          ";
    file << endl;
  }
}

friend ostream& operator<<(ostream& os, const complex& a)
 {
   os<<a.re<<" " << a.im<<"               ";
   return os;
 }

ostream& operator<<(ostream& os, cvect& a)
  {
    for(int i=0; i<N;i++)
      os<<a[i]<<"   ";
    return os;
  }
ostream& operator<<(ostream& os, cmat& a)
  {
    for(int i=0; i<N;i++)
     for(int j=0; j<N;j++)
        os<<a(i,j)<<"   ";
    return os;
  }

ここに N=2

2.17. メッシュの適応制御

adaptmeshの予約語を使うと、最初の計算結果に基づいて、メッシュの品質 を高めるための「適応制御」(mesh adaptation)を行うことができます。次の例を 見てください。

Navier stokes 方程式の後退ステップです。

nowait;
border(1,0,1,6,5)   begin x:=0;       y:=1-t  end;
border(2,0,1,15,2) begin x:=2*t;     y:=0     end;
border(2,0,1,10,2) begin x:=2;       y:=-t    end;
border(2,0,1,20,2) begin x:=2+3*t;   y:=-1    end;
border(2,0,1,35,2) begin x:=5+15*t; y:=-1     end;
border(4,0,1,10,2) begin x:=20;      y:=-1+2*t end;
border(5,0,1,35,5) begin x:=5+15*(1-t); y:=1  end;
border(5,0,1,40,5) begin x:=5*(1-t);y:=1      end;


buildmesh(1500);
nu := 0.001; dt := 0.4;


/* initial values */
u = y*(1-y)*one(y>0);
v = 0;
p = 2*nu*x*one(y>0);
un = u; vn = v;
i:=1; j:=2; k:=3;count:=1;
iter(500)
begin
 if (count % 100 == 0) then { /* adapt every 100 iterations */
   adaptmesh(u,p);
   i:=1; j:=2;k:=3;
   wait; /* wait after each displayed solution */
   plot(u);
   plot(v);
   plot(p);
   nowait; /* do not wait anymore */
 };
 /* increment the counter */
 count:=count+1;


 f=convect(un,u,v,dt);   g=convect(vn,u,v,dt);


 solve(u,i){          /*Horizontal velocity*/
   onbdy(1) u = y*(1-y);
   onbdy(2,5) u = 0;
   onbdy(4)dnu(u)=0;
   pde(u) id(u)/dt-laplace(u)*nu = f/dt -dx(p);
 };


 solve(v,j){         /* Vertical velocity */
   onbdy(1,2,4,5) v = 0;
   pde(v) id(v)/dt-laplace(v)*nu = g/dt -dy(p);
 };


 solve(p,k) {       /*   Pressure */
 onbdy(1,2,5) dnu(p) = 0;
 onbdy(4) p=0;
 pde(p)   -laplace(p)= -(dx(f) + dy(g))/dt;
 };
 un = u;   vn = v; i:=-1; j:=-2; k:=-3;
end ;
wait;
plot(u);plot(v);plot(p);

500 回の繰り返し中、100 回毎に適応制御を行っています。

  w=0;
  plot(w);

borderに新しい引数が追加されていることに注意してください。この最後 の引数(オプション)はこの境界の起点と見なされます。メッシュの適応制御を行う ときはこれを指定してください。さもないと、未定義の辺や頂点ができることがあ ります。

2.18. インタプリータの構造

インタプリータが使うtokenは lexical.h で定義され、字句解析は lexical.cxx の netsym 関数で行われます。

構文解析は syntaxic.cxx の instruction 関数が行いますが、構文は LL(1) で 設計されています。

構文解析プログラムは入力プログラムを plante 関数により構文解析木に変換し、 最後に eval 関数を使って、その root から実行を始めます。

FreeFEMに新しい命令を追加する場合は、次のルールを守る必要があります。

• 構文が LL(1) になっていて、既存の命令と矛盾しないこと
• 予約語を installe (lexical.cxx) に追加し、symbols (lexical.h) を更新して、 instruction に構文解析に対応する C の関数を追加し、eval の switch に追加する。

構文解析木のノードで使われている構造体は、下記のような極めて簡単なものです。

  typedef struct noeud {
    Symbol symb;
    creal value;
    ident *name;
    long junk;
    char *path;
    struct noeud *l1, *l2, *l3, *l4;
  } noeud, *arbre;

3. 例

3.1 三角分解(Triangulations examples)

o リング

  twopi := 2*pi;
  border(1,0,twopi,60) { x := cos(t); y := sin(t) };
  border(2,0,twopi,20) { x := 0.25*cos(-t); y := 0.25*sin(-t) };
  buildmesh(400);

o 長方形 [(0,0),(0,2),(2,1),(0,1)]

  border(1,0,2,20) { x:= t; y:= 0 };
  border(1,0,1,10) { x:= 2; y:= t };
  border(1,0,2,20) { x:= 2-t; y:= 1 };
  border(1,0,1,10) { x:= 0; y:= 1-t };
  buildmesh(400);

o 境界番号を細かく割り当てた正方形

  border(1,0,4,41)
  {
    if(t<=1) then { x:=t; y:=0 };
    if((t>1)and(t<2)) then { x:=1; y:=t-1; ib:= 2 };
    if((t>=2)and(t<=3)) then { x:=3-t; y:=1; ib:= 3 };
    if(t>3) then { x:=0; y:=4-t; ib:= 4 }
  };
  buildmesh(400);

o 複数の領域を持つ円

  border(1,0,2,17) { x:= cos(pi*t); y:= sin(pi*t)};
  border(0,-1,1,7) { x:= t; y:=0; };
  border(0,0,1,4) { x:=0;y:=t };
  buildmesh(300);
  /* observe the value of "region" by using "show triangle numbers" */

3.2. スカラ変数

o キャパシタ(電気コンデンサ)

  /* a circle of radius 5 centered at (0,0) */
  border(1,0,2*pi,60) { x:=5*cos(t); y:=5*sin(t) };
  /* The rectangle on the right */
  border(2,0,1,4) { x:=1+t; y:=3 };
  border(2,0,1,24) { x:=2; y:=3-6*t };
  border(2,0,1,4) { x:=2-t; y:=-3 };
  border(2,0,1,24) { x:=1; y:=-3+6*t };
  /* The rectangle on the left */
  border(3,0,1,4) { x:=-2+t; y:=3 };
  border(3,0,1,24) { x:=-1; y:=3-6*t };
  border(3,0,1,4) { x:=-1-t; y:=-3 };
  border(3,0,1,24) { x:=-2; y:=-3+6*t };
  buildmesh(800);

  /* Boundary conditions and PDE */
  solve(v) {
    onbdy(1) v = 0;
    onbdy(2) v = 1;
    onbdy(3) v = -1;
    pde(v) -laplace(v) = 0;
  };
  plot(v);

o 熱伝達と放射

  border(1,0,22,89)
  {
    if(t<=10)then { x:= t; y:=0 ; ib:=3 };
    if((t>10)and(t<11))then { x:=10; y:=t-10; ib:=2 };
    if((t>=11)and(t<=21))then { x:=21-t; y:=1; ib:=4 };
    if(t>21)then { x:=0; y:=22-t; ib:=1 };
  };
  buildmesh(800);

  changewait;
  t0 := 10;   t1 := 100;   te := 25; b=0.1; c = 5.0e-8;
  w = (b + 2*c * (te+546)*(te+273)*(te+273));
  solve(v,1) {
    onbdy(1) v=t0;
    onbdy(2) v = t1;
    onbdy(3) dnu(v)=0;
    onbdy(4) id(v) * w + dnu(v) = te * w;
    pde(v) -laplace(v) * y =0;
  };
  iter(10) {
    u=v;
    w = (b + c * (u+te + 546)*((u+273)*(u+273) + (te+273)*(te+273)));
    solve(v,-1) {
      onbdy(1) v=t0;
      onbdy(2) v = t1;
      onbdy(3) dnu(v)=0;
      onbdy(4) id(v)*w + dnu(v)= te * w;
      pde(v) -laplace(v) * y = 0
    };
  };
  plot(v);

対流による熱伝達の簡単なモデル

  (d/dn)(u) = -b*(u - te)

  (d/dn)(u) - c*((u + 273)^4 - (te + 273)^4)) = 0

  (d/dn)(u) + b*(u-te) - c*((u+273)^4 - (te+273)^4)) = 0

  a^4 - b^4 = (a - b)*(a + b)*(a^2 + b^2)

  (d/dn)(U) + b*(U-te) - c*(U-te)*(u+te+546)*((u+273)^2 + (te+273)^2)) = 0

  v = u - te

o 熱伝導: 一様でない材料

  h0:=1.0; h1:=2.0;
  border(1,0,22,89) {
    if(t<10)then { x:=t; y:=0; ib:=3 };
    if((t>=10)and(t<=11))then { x:=10; y:=h1*(t-10); ib:=2 };
    if((t>11)and(t<21))then { x:=21-t; y:=h1; ib:=4 };
    if(t>=21)then { x:=0; y:=h1*(22-t)i; ib:=1 };
  };
  border(0,0,10,41) { x:=t; y:=h0 };
  buildmesh(800);

  t0=10; t1=100; alpha:=1;
  kappa=0.01 + one(y>h0);
  solve(v) {
    onbdy(1) v=t0;
    onbdy(2) v=t1;
    onbdy(4) dnu(v)=0;
    onbdy(3) dnu(v)=0;
    pde(v) -laplace(v)*kappa*y + id(v)*kappa*y = 0;
  };
  plot(v);

o 圧縮性ポテンシャル流れ

  changewait;/* gamma = 1.4, outer circle radius is 5 */
  mach1 := 1/sqrt(6); machinfty = 0.85*mach1;
  rhoinfty=sqrt((1-machinfty^2)^5);
  solve(phi) begin
    onbdy(1) dnu(phi) = rhoinfty*machinfty*x/5; onbdy(2) dnu(phi) = 0;
    pde(phi) id(phi)*0.0001-laplace(phi) = 0;
  end;
  u1 = dx(phi); u2 = dy(phi); rho=sqrt((1-(u1^2+ u2^2))^5); plot(phi);
  iter(5)
  begin
    solve(phi) begin
      onbdy(1) dnu(phi)=rhoinfty*machinfty*x/5;
      onbdy(2) dnu(phi) = 0;
      pde(phi) id(phi)*0.0001-laplace(phi)*rhop=0;
    end;
    u1 = dx(phi); u2 = dy(phi); rho=sqrt((1-(u1^2+ u2^2))^5);
    rhop = convect(rho,u1,u2,0.1); plot(rho)
  end;
  plot(sqrt((u1^2+u2^2))/mach1);

o ナビエ・ストークス方程式

  /* Poor but better than none algorithm*/
  changewait;
  border(1,0,1,6)  { x:=0;          y:=1-t    };
  border(2,0,1,15) { x:=2*t;        y:=0      };
  border(2,0,1,10) { x:=2;          y:=-t     };
  border(2,0,1,20) { x:=2+3*t;      y:=-1     };
  border(2,0,1,35) { x:=5+15*t;     y:=-1     };
  border(3,0,1,10) { x:=20;         y:=-1+2*t };
  border(4,0,1,35) { x:=5+15*(1-t); y:=1      };
  border(4,0,1,40) { x:=5*(1-t);    y:=1      };
  buildmesh(900);

  nu = 0.002; dt := 0.4;

  /* initial pressure */
  solve(p,1) {
    onbdy(1)dnu(p) =-2*nu;
    onbdy(3) p=0; onbdy(2,4) dnu(p) = 0;
    pde(p) - laplace(p)= 0;
  };
  /* initial horizontal velocity */
  solve(u,2) {
    onbdy(1) u = y*(1-y);
    onbdy(3) dnu(u) = 0; onbdy(2,4) u = 0;
    pde(u) id(u)/dt-laplace(u)*nu = -min(y*y-y,0)/dt;
  };
  /* initial vertical velocity */
  solve(v,3){
    onbdy(1,3)v = 0; onbdy(2,4) v = 0;
    pde(v) id(v)/dt-laplace(v)*nu =0;
  };
  un = u; vn = v;
  iter(80)
  {
    f=convect(un,u,v,dt); g=convect(vn,u,v,dt);
    /*Horizontal velocity*/
    solve(u,-2) {
      onbdy(1) u = y*(1-y); onbdy(2,4) u = 0;
      onbdy(3)dnu(u)=0;
      pde(u) id(u)/dt-laplace(u)*nu = f/dt -dx(p);
    };
    plot(u);
    /* Vertical velocity */
    solve(v,-3) {
      onbdy(1,2,3,4) v = 0;
      pde(v) id(v)/dt-laplace(v)*nu = g/dt -dy(p);
    };
    /* Pressure */
    solve(p,-1) {
      onbdy(1)dnu(p) =-2*nu;
      onbdy(3) p=0; onbdy(2,4) dnu(p) = 0;
      pde(p) -laplace(p)= -(dx(f) + dy(g))/dt;
    };
    un = u; vn = v;
  };
  save('u.dta',u); save('v.dta',v); save('p.dta',p); plot(u);

3.3. 複素変数

  complex; nowait;
  border(1,0,1,10) { x:=t; y:=0; };
  border(1,0,1,10) { x:=1; y:=t; };
  border(2,0,1,10) { x:=1-t; y:=1; };
  border(1,0,1,10) { x:=0; y:=1-t; };
  buildmesh(200);

  solve(u)   /* observe than Re(u) = Im(u) */
  {
    onbdy(1,2) u=0;
    pde(u) id(u)-laplace(u)=1+I ;
  };
  v=Im(u);
  plot(u);plot(v);plot(u-v);

3.4. 連立方程式

  /* This is a 2-system example for which the solution is know
  analytically, thus the precision of Gfem can be checked */
  nowait;
  ns:=40;
  border(1,0,2*pi,2*ns) { x:= 3*cos(t); y:= 2*sin(t); };
  border(2,0,2*pi,ns)    { x:= cos(-t); y:= sin(-t); };
  buildmesh(ns*ns);

  ue= sin(x+y);
  ve = ue;
  p = ue;
  nx = -x;
  ny =- y;
  dxue = cos(x+y);

  c = 0.2;
  a1 = y;
  a2 =x;
  nu = 1;
  nu11 = 1;
  nu22 = 2;
  nu21 =0.3;
  nu12 =0.4;
  b=1;

  dnuue=dxue*(nu*(nx+ny) +
  (nu11 + nu12)*nx + (nu21+ nu22)*ny);
  g = ue*c+dnuue;
  f = b*ue+dxue*(a1+a2) +ue*(2*nu+nu11+nu12+nu21+nu22);

  solve(u,v) {
    onbdy(1) u = p;
    onbdy(1) v = p;
    onbdy(2) id(u)*c/2 + id(v)*c/2 + dnu(u) = g;
    onbdy(2) id(v)*c + dnu(v) = g;
    pde(u) id(u)*b + dx(u)*a1 + dy(u)*a2
             -laplace(u)*nu - dxx(u)*nu11 -
             dxy(u)*nu12 - dyx(u)*nu21 - dyy(u)*nu22 =f;

    pde(v) id(v)*b/2+id(u)*b/2
             + dx(v)*a1 + dy(v)*a2 -laplace(v)*nu
             - dxx(v)*nu11 - dxy(v)*nu12
             - dyx(v)*nu21 - dyy(v)*nu22 = f;
  };

  plot(abs(u-ue) + abs(v-ve));

4. 注

このドキュメントはFreeFEMのオリジナル・ドキュメント FreeFEM Documentationを書き換えたものです。説明不足と思われる部分は わかる範囲で補足しましたが、省略した部分もあります。

4.1 注 1

  (d/dt)A - (d/dx)(n*(d/dx)(A)) + (d/dy)(n*(d/dy)(A)) = g*V

  (Am+1 - Am)/Dt - (n/2)*laplace(Am+1 + Am) = g*V

  w = (Am+1 + Am)/2

  (2/Dt)*(w - Am) - n*lalplace(w) = g*V
  Am+1 = 2*w - Am

  solve(w) {
    pde(w) id(w)*2/dt-laplace(w)*n = g*V + a*2/dt
  };
  a = 2*w - a

5. ローカルな機能追加

以下、FreeFEMのオリジナルにはなくて、ローカル版で私が自分で使うために 追加した機能です。

save

save() のファイル名として「-」を指定すると、stdout に書き出します。結果の 表示に便利です。

数学関数

C と同じ、atan2(y,x), fmod(x,y), floor(x) を追加しました。

while() 文

while(logical_expr) 文を追加しました。

for() 文

for(instruction; logical_expr; instruction) を追加しました。

++, --, +=, -=, *=, /=, &&, || 演算子

「++」「--」「+=」「-=」「*=」「/=」「&&」「||」の演算子を追加しました。

delay() 関数の追加

delay(n) で n マイクロ秒の遅延が得られます。連続表示のタイミング調整に使います。

border() の機能

border(n,...) の n に定数だけでなく、式も使えるようにしました。

要素分割の領域(region)表示

「-r」オプションを指定するか要素数が少ないとき、region 番号も表示するよう にしました。

要素分割の境界(boundary)番号表示

「-b」オプションを指定すると境界番号も表示するようにしました。

printf(fmt,p1,p2,p3)

C と同じ printf() 文を追加しました。パラメータは最大３つまでしか使えません。

onbdy(b)

onbdy(b) の b に式も使えるようにしました。

plotmesh

「Introduction to Scientific Computing」で使われている plotmesh を追加 しました。

また、オリジナル版のバグと思える、次の点を修正しました。

  1) intt[u] が intt(0)[u] だけでなく、全ての要素を積分するようにした
  2) save('file',f) の f が実数のとき複素数形式で書かないようにした
  3) iter() のパラメータは定数だけでなく、式も使えるようにした
  4) 「%」演算子が fmod() として機能するようにしました
  5) 「!=」演算子の token を認識するようにした
  6) exec() 関数が実際に使えるようにした

別のテキストを見ると、既に 2) は未公開版で修正されているようです。