66 lines
1.9 KiB
Matlab
66 lines
1.9 KiB
Matlab
clc
|
||
clear all
|
||
|
||
%参考四面体
|
||
%顶点1-[1,0,0] 2-[0,1,0] 3-[0,0,1] 4-[0,0,1]
|
||
|
||
% 算例配置（改波长/材料请编辑 case_config.m）
|
||
[phy, cfg] = case_config();
|
||
|
||
%读取网格（与 C++ SBCmesh.dat 一致；无 .mat 时读 .dat）
|
||
rootDir = fileparts(mfilename('fullpath'));
|
||
matFile = fullfile(rootDir, cfg.meshMatFile);
|
||
datFile = fullfile(rootDir, cfg.meshDatFile);
|
||
if isfile(matFile)
|
||
load(matFile, 'mesh');
|
||
elseif isfile(datFile)
|
||
mesh = load_SBCmesh_dat(datFile);
|
||
else
|
||
altDat = fullfile(rootDir, '..', '..', '..', '3D opticsfem-master', 'SBCmesh.dat');
|
||
if isfile(altDat)
|
||
mesh = load_SBCmesh_dat(altDat);
|
||
else
|
||
error('未找到 SBCmesh.mat 或 SBCmesh.dat，请先运行 getMesh.m 或复制 C++ 网格。');
|
||
end
|
||
end
|
||
|
||
%矩阵组装
|
||
solver.Ai=[];solver.Aj=[];solver.Av=[];
|
||
solver.Bi=[];solver.Bj=[];solver.Bv=[];
|
||
solver=assembly_equ(phy,mesh,solver);
|
||
solver.A=sparse(solver.Ai,solver.Aj,solver.Av);
|
||
solver.B=sparse(solver.Bi,solver.Bj,solver.Bv);
|
||
|
||
% 导出 A/B 矩阵（与 C++ FemType=5 OutFile 格式一致）
|
||
export_ab_matrix(solver);
|
||
|
||
%求解设置
|
||
targetff0=phy.ff0;
|
||
targetlam0=phy.c_const/targetff0;
|
||
targetk0=2*pi/targetlam0;
|
||
sigma=targetk0*targetk0;
|
||
numberSolve=cfg.NbrMode;
|
||
%求解
|
||
[x,GAMA_sq]=eigs(solver.A,solver.B,numberSolve,sigma);
|
||
solverk0=diag(sqrt(GAMA_sq));
|
||
solverlam0=2*pi./solverk0;
|
||
solverff0=phy.c_const./solverlam0
|
||
|
||
|
||
%计算电场
|
||
Ex=zeros(mesh.NbrVertex,numberSolve);
|
||
Ey=zeros(mesh.NbrVertex,numberSolve);
|
||
Ez=zeros(mesh.NbrVertex,numberSolve);
|
||
normE=zeros(mesh.NbrVertex,numberSolve);
|
||
for i=1:numberSolve
|
||
xx=x(:,i);
|
||
[Ex(:,i),Ey(:,i),Ez(:,i),normE(:,i)]=get_ele(mesh,xx);
|
||
end
|
||
|
||
% 导出 freq / Ex / Ey / Ez / normE（与 C++ Post_3D_EigenFreq 格式一致）
|
||
export_eigen_fields(solverff0, Ex, Ey, Ez, fullfile(fileparts(mfilename('fullpath')), 'OutFile'));
|
||
|
||
%绘图
|
||
triIndex=findTri(3,mesh);
|
||
tri=mesh.Tri(triIndex,:);
|
||
plotE(mesh.Vertex,tri,normE(:,3),0); |