openEMSのメッシュを切る関数自作した

はじめに

openEMSはフリーでオープンソースのFDTD（有限差分時間領域法）の電磁界解析ソフトで、僕は結構気に入って使っているのだが、日本ではあまり流行っていない。もっと広まってもいいと思うんだけどなあ。

openEMSは、MatlabまたはGNU Octaveで3Dモデルを作成する。他の有料の電磁界解析ソフト（HFSSとかCOMSOLとか）と違って、メッシュも自分で切らなければいけないところが難しい。一応、メッシュをいい感じに切ってくれるSmoothMesh関数などが用意されているが、これがなかなか使いにくい。処理が早く終わる時もあれば条件によっては延々と終わらない時もあって、結構癖のある関数だ。

そこで、簡単にopenEMSのメッシュを分割する関数を自作してみた。ソースコードは一番下に載せる。

メッシュ分割の考え方

おおざっぱなイメージを持ってもらうために下の例を見てほしい。

メッシュ分割の例

FR4の基板上にPEC（perfect electric conductor）で作ったマイクロストリップ線路と中央にDUTがあるモデル。メタルのエッジには電界が集中しやすいのでメッシュは細かく切りたい。一方で、オブジェクトから離れた電界が弱い部分はメッシュは粗くても問題ない。エッジ付近は細かくして、徐々に間隔を広げていくようなメッシュの切り方が望ましい。全部細かく分割してしまうのは、無駄に計算時間が長くなってしまい効率的でない。

上のようなメッシュ分割法をもう少し具体的に表すと、

• オブジェクトのエッジ周辺は最小メッシュ分解能（min_res）で分割する
• 一定の比率（ratio）で徐々にメッシュの間隔を大きくする
• 最大メッシュ分解能（max_res）に達したら、等間隔に分割する

ここで、以下のパラメータを定義しておく。max_res：最大のメッシュ分解能。min_res：最小のメッシュ分解能。ratio：隣り合うメッシュの比の上限値。

目安として、max_resは波長λの1/30、min_resはmax_res/5、ratioは1.3くらいに設定するとよいと思う。

では、実際にx1からx2の区間（x2>x1）を上のルールに従って分割する方法を考えてみよう。ここでは、x=x1にオブジェクトのエッジがあると仮定して、x=x1近傍は細かく、x=x2に向かうほど間隔が広がっていくようにしたい。

といっても、僕は数学的な難しいことは苦手なので、小学生でも思いつきそうな方法でやってみる。まず、長さが (x2-x1) のテープを用意する。

テープの左端からmin_resの長さ分だけハサミで切り取る。

ratioの比率で徐々に切る長さを増やしていく。長さがmax_resに達したら、あとはmax_resで等間隔に切っていく。

最後に切れ端（あまり）が残る。切れ端の長さで場合分けする。

i) 切れ端の長さ ≥ min_res のとき
これは簡単。小さい順に並び替えるだけ（ソート）。並び替えた後の切れ端と隣接するメッシュとの比は必ずratio以下になるのでOK。

ii) 切れ端の長さ < min_res のとき
ii-a) 右二つの長さの合計 ≥ 2×min_res のとき
切れ端と右から2番目（つまり一番長いやつ）を一旦くっつけて、切れ端がmin_resになるように再分配する。あとはソートすればOK。

ii-b) 右二つの長さの合計 < 2×min_res のとき
この場合分けは、分割数が少ないときに発生する。この場合、切れ端をゴミ箱にぽいして、残りのテープの長さの合計が (x2-x1) になるようにスケーリングする。

以上のアルゴリズムでメッシュ分割できるようになった。あとは実装するだけ。

実装・関数の使い方

Matlab/GNU Octaveの関数として実装した。ソースコードは一番下に載せた。ここでは、関数の使い方を説明する。

GradedMeshLines()

GradedMeshLines()は1次元のメッシュ分割をする関数。

function lines = GradedMeshLines(lines, max_res, min_res, ratio)

引数に、もととなるグリッド（両端の座標、オブジェクトのエッジの座標の配列）と、max_res、min_res、ratioを渡す。min_resとratioはオプションでそれぞれmax_res/5、1.3をデフォルト値としている。

使い方の例↓（そのまま実行可能）

lines = [-10 0 10]; % もととなるグリッド
lines = GradedMeshLines(lines, 1); % メッシュ分割

% プロット
figure('Position', [100 100 600 100]);
for x = lines
    plot([x x], [0 1], 'k-');
    hold on;
end

これの実行結果が下。

GradedMeshLines()の使い方

多くの場合、シミュレーションモデルの両端は吸収境界に割り当てるので、両端はグレーディングせず、内側のラインをオブジェクトのエッジとみなしてメッシュ分割するようにしてある。

GradedMesh()

GradedMesh()は3次元のメッシュ分割をする関数。

function mesh = GradedMesh(mesh, max_res, min_res, ratio)

openEMSのMatlab/Octave関数の中にDetectEdges()という便利な関数があるので、最初にDetectEdges()でCSXのオブジェクトのエッジを検出してから、GradedMesh()を適用すれば簡単にメッシュが作成される。

使い方の例↓（openEMSのCSXオブジェクト作成後）

resolution = 1; % mm
mesh = DetectEdges(CSX); % CSXオブジェクトのエッジを検出
mesh = GradedMesh(mesh, resolution); % メッシュを分割
CSX = DefineRectGrid(CSX, unit, mesh); % 直交座標系のグリッドを設定

ソースコード

ファイル名：GradedMeshLines.m

function lines = GradedMeshLines(lines, max_res, min_res, ratio)
% function lines = GradedMeshLines(lines, max_res, min_res, ratio)
%
% input parameters:
%   lines                   lines to create a graded mesh in between.
%   max_res                 maximum resolution.
%   min_res (optional)      minimum resolution. default is max_res/5.
%   ratio (optional)        ratio of adjacent grids. default is 1.3.

    if ~exist('min_res', 'var')
        min_res = max_res/5;
    end
    if ~exist('ratio', 'var')
        ratio = 1.3;
    end

    lines = uniquetol(lines);

    if length(lines) < 3
        n = 1+floor(abs(lines(end)-lines(1)));
        lines = linspace(lines(1), lines(end), n);
    else
        lines_copy = lines;
        lines = [lines singleGradedMesh(lines_copy(2), lines_copy(1), max_res, min_res, ratio)];
        lines = [lines singleGradedMesh(lines_copy(end-1), lines_copy(end), max_res, min_res, ratio)];
        for n = 2:length(lines_copy)-2
            lines = [lines doubleGradedMesh(lines_copy(n), lines_copy(n+1), max_res, min_res, ratio)];
        end
        lines = sort(lines);
    end

end

function lines = singleGradedMesh(x1, x2, max_res, min_res, ratio)
    delta_total = abs(x1-x2);
    delta = [min_res];

    while sum(delta) < delta_total
        delta = [delta min(ratio*delta(end), max_res)];
    end
    delta(end) = delta_total-sum(delta(1:end-1));

    if delta(end) < min_res && length(delta) > 1
        if sum(delta(end-1:end)) > 2*min_res
            delta(end-1) = sum(delta(end-1:end))-min_res;
            delta(end) = min_res;
        else
            scaling = delta_total/sum(delta(1:end-1));
            delta = scaling*delta(1:end-1);
        end
    end

    delta = sort(delta);
    lines = [];
    for n = 1:length(delta)-1
        lines = [lines x1+sign(x2-x1)*sum(delta(1:n))];
    end
end

function lines = doubleGradedMesh(x1, x2, max_res, min_res, ratio)
    delta_total = abs(x1-x2);
    delta = [min_res];

    while sum(delta) < delta_total
        if mod(length(delta), 2) == 0
            delta = [delta min(ratio*delta(end), max_res)];
        else
            delta = [delta delta(end)];
        end
    end
    delta(end) = delta_total-sum(delta(1:end-1));
    
    if delta(end) < min_res && length(delta) > 1
        if sum(delta(end-1:end)) > 2*min_res
            delta(end-1) = sum(delta(end-1:end))-min_res;
            delta(end) = min_res;
        else
            scaling = delta_total/sum(delta(1:end-1));
            delta = scaling*delta(1:end-1);
        end
    end

    delta = sort(delta);
    delta = [delta(1:2:end) flip(delta(2:2:end))];

    lines = [];
    for n = 1:length(delta)-1
        lines = [lines x1+sign(x2-x1)*sum(delta(1:n))];
    end
end

ファイル名：GradedMesh.m

function mesh = GradedMesh(mesh, max_res, min_res, ratio)
    if ~exist('min_res', 'var')
        min_res = max_res/5;
    end
    if ~exist('ratio', 'var')
        ratio = 1.3;
    end
    
    mesh.x = GradedMeshLines(mesh.x, max_res, min_res, ratio);
    mesh.y = GradedMeshLines(mesh.y, max_res, min_res, ratio);
    mesh.z = GradedMeshLines(mesh.z, max_res, min_res, ratio);

end