115 lines
3.0 KiB
Matlab
115 lines
3.0 KiB
Matlab
function [Solver,Mesh]=FemMatrixAssembly(Physic,Material,Mesh,...
|
|
Port1Mesh,Port1Solver,Port2Mesh,Port2Solver)
|
|
|
|
|
|
%process maxwell's equations
|
|
tic
|
|
epsilonr=Material.epsilonr-1i*Material.sigma*Material.temp;
|
|
[A,Mesh.Edges,Mesh.EdgesOfElements,Dof]=PhysicMatrixAssembly(Mesh.Nodes,Mesh.Elements,...
|
|
Mesh.Domains,Physic.lam0,epsilonr,Material.mur);
|
|
b=complex(zeros(Dof,1));
|
|
disp('The time taken of matrix assmebly for maxwell equations');
|
|
toc
|
|
|
|
%process boundary mesh
|
|
tic
|
|
el2no=Mesh.Elements';
|
|
f1=el2no([1 1 1 2],:);
|
|
f2=el2no([2 2 3 3],:);
|
|
f3=el2no([3 4 4 4],:);
|
|
f_array=[f1(:) f2(:) f3(:)];
|
|
[Faces,~,FacesOfElements]=unique(f_array,'rows');
|
|
Mesh.AllFaces=Faces;
|
|
Mesh.FacesOfElements=FacesOfElements;
|
|
addDof=0;
|
|
|
|
%process port1 boundary condition
|
|
Port1Mesh.FacesConn=GetBFacesConn(Port1Mesh.OldFaces,Faces,FacesOfElements);
|
|
[A,b,addDof]=PortBCMatrixAssembly(A,b,Mesh.Nodes,Mesh.Elements,Mesh.EdgesOfElements,...
|
|
Mesh.Domains,Material.mur,Port1Mesh,Port1Solver,Dof,addDof);
|
|
|
|
%process port2 boundary condition
|
|
Port2Mesh.FacesConn=GetBFacesConn(Port2Mesh.OldFaces,Faces,FacesOfElements);
|
|
[A,b,addDof]=PortBCMatrixAssembly(A,b,Mesh.Nodes,Mesh.Elements,Mesh.EdgesOfElements,...
|
|
Mesh.Domains,Material.mur,Port2Mesh,Port2Solver,Dof,addDof);
|
|
addDof(1)=[];
|
|
disp('The time taken of processing port boundary condition');
|
|
toc
|
|
|
|
%process PEC boundary condition
|
|
[PECFaces,~]=GetBFaces(Mesh.Faces,Mesh.FacesIndex,Physic.PEC);
|
|
fl2no=PECFaces';
|
|
n1=fl2no([1 1 2],:);
|
|
n2=fl2no([2 3 3],:);
|
|
fl_fd2no_array=[n1(:) n2(:)];
|
|
PECEdges=unique(fl_fd2no_array,'rows');
|
|
PECInd=find(ismember(Mesh.Edges,PECEdges,'rows'));
|
|
Free_ind=1:(Dof+sum(addDof));
|
|
Free_ind(PECInd)=[];
|
|
A=A(Free_ind,Free_ind);
|
|
b=b(Free_ind);
|
|
%Matrix reordering
|
|
r=symrcm(A);
|
|
A=A(r,r);
|
|
b=b(r);
|
|
|
|
%solving
|
|
tic
|
|
[P,R,C] = equilibrate(A);
|
|
newA=R*P*A*C;
|
|
newb=R*P*b;
|
|
|
|
xx=newA\newb;
|
|
xx=C*xx;
|
|
|
|
% [Ai,Aj,As]=find(A);
|
|
% dlmwrite('Ai.txt',Ai,'delimiter',' ','precision',20);
|
|
% dlmwrite('Aj.txt',Aj,'delimiter',' ','precision',20);
|
|
% dlmwrite('As.txt',As,'delimiter',' ','precision',20);
|
|
% dlmwrite('b.txt',b,'delimiter',' ','precision',20);
|
|
|
|
% fileID = fopen('x_real_0.txt');
|
|
% C1 = textscan(fileID,'%f');
|
|
% fclose(fileID);
|
|
% xx=C1{:};
|
|
% fileID = fopen('x_complex_0.txt');
|
|
% C1 = textscan(fileID,'%f');
|
|
% fclose(fileID);
|
|
% xx=xx+1i*C1{:};
|
|
% xx=C*xx;
|
|
|
|
%processing solution
|
|
q=zeros(length(r),1);
|
|
for i=1:length(r)
|
|
q(r(i))=i;
|
|
end
|
|
xx=xx(q);
|
|
%processing PEC
|
|
for i=1:length(PECInd)-1
|
|
n=PECInd(i);
|
|
xx=[xx(1:n-1,:);zeros(1,1);xx(n:end,:)];
|
|
end
|
|
i=length(PECInd);
|
|
n=PECInd(i);
|
|
if n<=length(xx)
|
|
xx=[xx(1:n-1,:);zeros(1,1);xx(n:end,:)];
|
|
else
|
|
xx=[xx(1:n-1,:);zeros(1,1);];
|
|
end
|
|
%compute the electric field
|
|
% Solver.condA=condest(newA);
|
|
Solver.Et=xx(1:Dof);
|
|
Solver.Dof=Dof;
|
|
Solver.e1=xx(Dof+1);
|
|
Solver.e2=xx(Dof+2);
|
|
disp('The time taken of solving');
|
|
toc
|
|
|
|
%Reconstructing electric field Ex,Ey,Ez and normE
|
|
tic
|
|
[Solver.Ex,Solver.Ey,Solver.Ez,Solver.normE]=GetExEyEz(1,...
|
|
Mesh.Nodes,Mesh.Elements,Mesh.EdgesOfElements,Solver.Et);
|
|
disp('The time taken of post-processing');
|
|
toc
|
|
|
|
end |