55 lines
1.8 KiB
Matlab
55 lines
1.8 KiB
Matlab
function solver = assembly_bele(phy, mesh, solver)
|
||
%ASSEMBLY_BELE BELE 背景场源项：方法 A 直接体积分
|
||
%
|
||
% b -= ∫ J_b·N dV，J_b = curl curl E_b − k₀²ε E_b （弱形式 Tt − St）
|
||
% 不拆 Tt−St，不在 SBC 外表面做 IBP / Be(E_b)
|
||
|
||
k0 = 2 * pi / phy.lda0;
|
||
|
||
[u, v, w, weight, nbrGP] = getGaussPoints(1);
|
||
|
||
for n = 1:mesh.NbrTet
|
||
domain = mesh.DomainOfTet(n);
|
||
if ~ismember(domain, phy.beleDomains)
|
||
continue;
|
||
end
|
||
|
||
x = mesh.Vertex(mesh.Tet(n, :), 1);
|
||
y = mesh.Vertex(mesh.Tet(n, :), 2);
|
||
z = mesh.Vertex(mesh.Tet(n, :), 3);
|
||
|
||
Jac = zeros(3, 3);
|
||
Jac(1, 1) = x(1) - x(4); Jac(1, 2) = y(1) - y(4); Jac(1, 3) = z(1) - z(4);
|
||
Jac(2, 1) = x(2) - x(4); Jac(2, 2) = y(2) - y(4); Jac(2, 3) = z(2) - z(4);
|
||
Jac(3, 1) = x(3) - x(4); Jac(3, 2) = y(3) - y(4); Jac(3, 3) = z(3) - z(4);
|
||
DetJac = abs(det(Jac));
|
||
|
||
E = zeros(3, 6, nbrGP);
|
||
for gp = 1:nbrGP
|
||
for j = 1:6
|
||
E(:, j, gp) = getBF(1, j, u(gp), v(gp), w(gp));
|
||
E(:, j, gp) = Jac \ E(:, j, gp);
|
||
end
|
||
end
|
||
|
||
eps = phy.eps(domain);
|
||
|
||
for i = 1:6
|
||
contrib = 0;
|
||
for gp = 1:nbrGP
|
||
px = x(4) + Jac(1, 1) * u(gp) + Jac(2, 1) * v(gp) + Jac(3, 1) * w(gp);
|
||
py = y(4) + Jac(1, 2) * u(gp) + Jac(2, 2) * v(gp) + Jac(3, 2) * w(gp);
|
||
pz = z(4) + Jac(1, 3) * u(gp) + Jac(2, 3) * v(gp) + Jac(3, 3) * w(gp);
|
||
bE = phy.bele.Eb(px, py, pz);
|
||
bCurlCurlE = phy.bele.curlcurlEb(px, py, pz);
|
||
Jb = bCurlCurlE - k0^2 * eps * bE;
|
||
contrib = contrib + weight(gp) * DetJac * dot(E(:, i, gp), Jb);
|
||
end
|
||
edgeId = mesh.EdgeOfTet(n, i);
|
||
% Weak-form BELE: b += Tt - St (same as 2D/C++ after IPP with A = St - Tt)
|
||
solver.b(edgeId) = solver.b(edgeId) - contrib;
|
||
end
|
||
end
|
||
|
||
end
|