Sample Scripts from GB Books
GB009:  3次元の電場

FlexPDEで3次元の場を扱う場合には、x-y平面に定義されたbase planeをz軸方向に押出す（extrusion）という考え方が基本になります。

右の図は電場の様子を調べようとする
3次元空間中の立方体ドメインを示しています。Base planeであるx-y平面上には (0, 0.5) の位置に +q の電荷（赤の点）が、(0, -0.5) の位置に -q の電荷（青の点）が置かれているものとし（q = 1e-10クーロン）、このときの周囲の電場をFlexPDEで計算します。この場合、静電ポテンシャル U を与える数式は明示的な形で与えることができるため（数式 (1) 参照）、偏微分方程式を解くというプロセスは伴いません。

点 (x₀, y₀, z₀ ) に電荷 q が置かれているときの静電ポテンシャル U を与える数式：

まずタイトルを設定します。
  TITLE
    'Two Point Charges of Different Signs'    { 3dfields1a.pde }

次に座標系が3次元直交座標系であることを明示します（デフォルトはCartesian2）。
  COORDINATES
    Cartesian3

関連する数式を定義します。静電ポテンシャル U を規定する数式もEQUATIONSセクションではなくDEFINITIONSセクション内で定義することになる点に注意してください。
  DEFINITIONS
    L = 1.0  d0 = 0.5  q = 1e-10
    eps0 = 8.85e-12
    c = 1/(4*PI*eps0)
    U = -q*c/sqrt(x^2+(y+d0)^2+z^2) + q*c/sqrt(x^2+(y-d0)^2+z^2)
                                          { Electrostatic potential }
    Ex = -dx(U)  Ey = -dy(U)  Ez = -dz(U) { Field components }
    E = -grad(U)  Em = magnitude(E)
    div_xy = dx(Ex) + dy(Ey)

z軸方向へのextrusionを定義します。下から順に境界面（Surface）とそれらにはさまれた層（Layer）を規定して行きます。今回の例では空間の属性がz軸方向で変化するわけではないのでレイヤの定義は必要ありません。
  EXTRUSION    { Extrude a cube through the (x, y) plane }
    Surface 'bottom' z = -L  { Limiting surfaces }
    Surface 'top' z = L

次にbase plane上で1辺が2の正方形ドメインを定義します。今回の例ではbase plane上の2次元ドメインが一様な形で -1 <= z <= 1 の範囲にextrudeされて行く形となります。
  BOUNDARIES
    Region 1
      Start(-L, -L)          { Trace outer boundary of base plane }
      Line to (L, -L) to (L, L) to (-L, L) to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y, z)
    Contour(U) painted on z = 0
    Contour(U) painted on z = 0.3
    Contour(U) painted on z = 1.0
    Vector(E) norm on z = 0
    Vector(E) norm on z = 1.0
    Elevation(Ex, Ey, Ez) from (0, -L, 1) to (0, L, 1)
    Contour(Ez) painted on z = 1.0 Report(Val(Ez, 0, 0.84, 1))
    Contour(div_xy) on z = 1.0
    Contour(div(E)) on z = 1.0

  END

(1) Grid(x, y, z)
FlexPDEによって生成された3次元のメッシュ構成を示しています。

(2) Contour(U) painted on z = 0
平面 z = 0 上におけるポテンシャル U の等高線図です。(0, 0.5, 0), (0, -0.5, 0) において U の理論値は無限大となるため、等高線図といっても曲線が極の周囲に局在する結果となっています。

(3) Contour(U) painted on z = 0.3
z = 0.3 という平面上には極は存在せず、U の値の範囲も限られたものとなるため、等高線図の形状も穏当なものとなっています。

(4) Contour(U) painted on z = 1.0
プロット(3)に比べると電荷の位置からさらに離れるため、曲面の形状はよりなだらかなものとなっています。

(5) Vector(E) norm on z = 0
平面 z = 0 上における電場 E のベクトルプロットを示したものです。normを指定しているため、場の強さはカラーによって表現されるのですが、ここでも極の存在によってモノトーンなプロットとなってしまっています。

(6) Vector(E) norm on z = 1.0
平面 z = 1 上での電場 E に関するベクトルプロットです。(0, 0.84, 1) の近傍に湧き出し口が、(0, -0.84, 1) の近傍に吸い込み口があるように見えます。0.84という値については次のelevationプロットを参照してください。

(7) Elevation(Ex, Ey, Ez) from (0, -L, 1) to (0, L, 1)
これは直線 (0, y, 1) 上におけるEx, Ey, Ezの変化をプロットしたものです。Exはどの場所でも 0 ですが、Eyは y = 0.84, -0.84 の位置で 0 となっています。つまり (0, 0.84, 1), (0, -0.84, 1) の2点では Ex = Ey = 0 となっているわけです。ただしこれらの点においてEzは 0 ではない点に注意してください。

(8) Contour(Ez) painted on z = 1.0 Report(Val(Ez, 0, 0.84, 1))
平面 z = 1上でのEzに関する等高線図です。グラフよりもReport文の出力に注意してください。(0, 0.84, 1) の点におけるEzの値が0.57であることがわかります。

(9) Contour(div_xy) on z = 1.0
平面 z = 1 上における

の値の等高線図です。これは 0 にはなりません。

(10) Contour(div(E)) on z = 1.0
一方、平面 z = 1 上における

の値は至るところ 0 となります。