音響数値解析のお勉強（有限要素法）

若手の音響研究者たちが作成しているOpenAcousticsというドキュメント＋ソースコード群があります。このサイトでは音響数値解析のソースコードだけでなく、充実した日本語の入門書が無料公開されています。だいぶ前に教えてもらっていたのですが、重い腰を上げて、空き時間に読み進めることにしました。

チュートリアルに使われているプログラミング言語はPythonです。Pythonはコードが理解しやすいので好きなのですが、より使い慣れているMatlabを使って自分用に書き直してみました。オリジナル版はGmshを使ってメッシュデータを読み込んでくるなど、より便利に使えるようになっています。下のコードではその部分は省略し、とりあえず動かしています。今度はドキュメントの数式を理解すべく、一行一行読んでいこうと思います。早くBEM、FDTD、CIPと進みたい。

（右のグラフは波の位置・周波数を変えて作ったもので、グラフの色もデフォルトから変更しています。）

function test_fem
c = 343.0;   % speed of sound

%% create triangular mesh grid of size N
N = 9;
[x,y] = meshgrid(linspace(0, 1, N), linspace(0, 1, N));
nodes = [x(:) y(:)];
elements = [];
for row = 1:length(x)-1;
  for col = 1:length(y)-1;
    n = (row-1)*length(x) + col;
    elements = [elements ; [n n+1 n+length(x)]];
    elements = [elements ; [n+1 n+length(x) n+length(x)+1]];
  end
end
numE = size(elements, 1);
numN = size(nodes, 1);

%% point source at node 1 with stationary state of 170Hz
q = zeros(numN, 1);
q(1) = 1;
f = 170;

%% compute and sum up element matrices
K = zeros(numN, numN);
M = zeros(numN, numN);
for n=1:numE
  [Kn,Mn] = triangularElement(nodes(elements(n,1),:), nodes(elements(n,2),:), nodes(elements(n,3),:));
  for d1=1:size(Kn, 1);
    for d2=1:size(Kn, 2);
      K(elements(n,d1), elements(n,d2)) = K(elements(n,d1), elements(n,d2)) + Kn(d1,d2);
      M(elements(n,d1), elements(n,d2)) = M(elements(n,d1), elements(n,d2)) + Mn(d1,d2) / c^2;
    end
  end
end

%% solve for phi
omega = 2 * pi * f;
phi = (K - (omega ^ 2) * M) \ q;

%% visualize
for n=1:numE
  patch(nodes(elements(n,:),1), nodes(elements(n,:),2), phi(elements(n,:)));
end
end


%% compute element matrix
function [K,M] = triangularElement(node0, node1, node2)
a = node1 - node0;
b = node2 - node1;
c = node0 - node2;

t0 = acos(-c * a' / (norm(c) * norm(a)));
t1 = acos(-a * b' / (norm(a) * norm(b)));
t2 = acos(-b * c' / (norm(b) * norm(c)));

s = (norm(a) + norm(b) + norm(c)) / 2;
S = sqrt(s * (s-norm(a)) * (s-norm(b)) * (s-norm(c)));

K = [
  cot(t1)+cot(t2) -cot(t2)         -cot(t1)
  -cot(t2)          cot(t0)+cot(t2) -cot(t0)
  -cot(t1)         -cot(t0)          cot(t0)+cot(t1)] / 2;

M = [
  2 1 1
  1 2 1
  1 1 2] * S / 12;
end