242 lines
7.6 KiB
Matlab
242 lines
7.6 KiB
Matlab
function [A,b,EdgesOfElements]=FemAssemble2(node,Elements,edge,edgeflag,domain,lam0,epsilonr,mur)
|
|
%圆形PML矩阵组装
|
|
|
|
%物理参数
|
|
k0=2*pi/lam0;
|
|
|
|
%高斯积分点
|
|
[xt,yt,pt,IntegralOrder]=GetGuassPoints(4);
|
|
% [lx,lp,IntegralEdgeOrder]=GetGuassLinePoints(4);
|
|
|
|
%网格处理
|
|
NbrElement=max(size(Elements));
|
|
el2no=Elements;
|
|
n1=el2no([1 1 2],:);
|
|
n2=el2no([2 3 3],:);
|
|
el_ed2no_array=[n1(:) n2(:)];
|
|
[Edges,~,EdgesOfElements]=unique(el_ed2no_array,'rows');
|
|
NbrEdges=max(size(Edges));
|
|
NbrEdge=max(size(edgeflag));
|
|
%边所在网格
|
|
ElementNum=zeros(NbrEdge*2,3);%1为所在网格 2为第几条边 2被边数
|
|
BNodes=edge';%对BNodes排序
|
|
for i=1:NbrEdge
|
|
if BNodes(i,1)>BNodes(i,2)
|
|
a= BNodes(i,1);
|
|
BNodes(i,1)=BNodes(i,2);
|
|
BNodes(i,2)=a;
|
|
end
|
|
end
|
|
ii=0;
|
|
for i=1:NbrEdge
|
|
for j=1:NbrElement*3
|
|
if (BNodes(i,1)==el_ed2no_array(j,1))&&(BNodes(i,2)==el_ed2no_array(j,2))
|
|
ii=ii+1;
|
|
ElementNum(ii,1)=fix((j-1)/3)+1;
|
|
ElementNum(ii,2)=j-(ElementNum(ii,1)-1)*3;
|
|
ElementNum(ii,3)=i;
|
|
end
|
|
end
|
|
end
|
|
|
|
%矩阵初始化
|
|
Dof=NbrEdges;
|
|
A=zeros(Dof,Dof);
|
|
b=zeros(Dof,1);
|
|
|
|
|
|
%循环网格 组装矩阵
|
|
for n=1:NbrElement
|
|
%节点坐标
|
|
x=zeros(3,1);y=zeros(3,1);l=zeros(3,1);
|
|
x(1)=node(Elements(1,n),1);y(1)=node(Elements(1,n),2);
|
|
x(2)=node(Elements(2,n),1);y(2)=node(Elements(2,n),2);
|
|
x(3)=node(Elements(3,n),1);y(3)=node(Elements(3,n),2);
|
|
l(1)=sqrt((x(1)-x(2))*(x(1)-x(2))+(y(1)-y(2))*(y(1)-y(2)));
|
|
l(2)=sqrt((x(1)-x(3))*(x(1)-x(3))+(y(1)-y(3))*(y(1)-y(3)));
|
|
l(3)=sqrt((x(3)-x(2))*(x(3)-x(2))+(y(3)-y(2))*(y(3)-y(2)));
|
|
%雅可比矩阵
|
|
Jac=zeros(3,3);
|
|
Jac(1,1)=x(2)-x(1);Jac(1,2)=y(2)-y(1);
|
|
Jac(2,1)=x(3)-x(1);Jac(2,2)=y(3)-y(1);Jac(3,3)=1;
|
|
InvJac=inv(Jac);
|
|
DetJac=abs(det(Jac));
|
|
DetJacx=(det(Jac));
|
|
TJac=Jac'/DetJacx;
|
|
%小矩阵
|
|
Et=zeros(3,3,IntegralOrder);curlEt=zeros(3,3,IntegralOrder);temp=zeros(3,1);
|
|
Ei=zeros(3,IntegralOrder);curlEi=zeros(3,IntegralOrder);
|
|
for i=1:IntegralOrder
|
|
for j=1:3
|
|
[temp(1),temp(2)]=BF_Et(j,xt(i),yt(i));temp(3)=0;
|
|
Et(j,:,i)=InvJac*temp*l(j);
|
|
temp(3)=BF_curlEt(j,xt(i),yt(i));temp(1)=0;temp(2)=0;
|
|
curlEt(j,:,i)=TJac*temp*l(j);
|
|
end
|
|
end
|
|
if domain(n)==5
|
|
for i=1:IntegralOrder
|
|
tempx=x(1)+Jac(1,1)*xt(i)+Jac(2,1)*yt(i);
|
|
tempy=y(1)+Jac(1,2)*xt(i)+Jac(2,2)*yt(i);
|
|
Ei(:,i)=Einc(tempx,tempy);
|
|
curlEi(:,i)=CurlEinc(tempx,tempy);
|
|
end
|
|
end
|
|
|
|
|
|
xxepsilonr=zeros(3,3,IntegralOrder);
|
|
xxinvXmur=zeros(3,3,IntegralOrder);
|
|
if domain(n)==5
|
|
xepsilonr=[epsilonr,0,0;0,epsilonr,0;0,0,epsilonr];
|
|
xmur=[mur,0,0;0,mur,0;0,0,mur];
|
|
invXmur=inv(xmur);
|
|
for i=1:IntegralOrder
|
|
xxepsilonr(:,:,i)=xepsilonr;
|
|
xxinvXmur(:,:,i)=invXmur;
|
|
end
|
|
else
|
|
for i=1:IntegralOrder
|
|
averX=x(1)+Jac(1,1)*xt(i)+Jac(2,1)*yt(i);
|
|
averY=y(1)+Jac(1,2)*xt(i)+Jac(2,2)*yt(i);
|
|
rho=sqrt(averX*averX+averY*averY);
|
|
rho0=1.5;
|
|
d=0.5;
|
|
R0=10;
|
|
sigma=((rho-rho0)/d)^2*R0/d/k0;
|
|
s1=1-1i*d/2/rho*sigma;
|
|
s2=1-1i*sigma;
|
|
aa=s1/s2;bb=s2/s1;cc=s1*s2;
|
|
Lambda=[(aa*averX*averX+bb*averY*averY)/rho/rho,((aa-bb)*averX*averY)/rho/rho,0;
|
|
((aa-bb)*averX*averY)/rho/rho,(bb*averX*averX+aa*averY*averY)/rho/rho,0;
|
|
0,0,cc];
|
|
xepsilonr=epsilonr*Lambda;
|
|
xmur=mur*Lambda;
|
|
invXmur=inv(xmur);
|
|
|
|
xxepsilonr(:,:,i)=xepsilonr;
|
|
xxinvXmur(:,:,i)=invXmur;
|
|
end
|
|
end
|
|
|
|
Sz=zeros(3,3);Tz=zeros(3,3);Bz=zeros(3,1);
|
|
for i=1:3
|
|
for j=1:3
|
|
for k=1:IntegralOrder
|
|
Sz(i,j)=Sz(i,j)+pt(k)*DetJac*dot(curlEt(i,:,k)',xxinvXmur(:,:,k)*curlEt(j,:,k)');
|
|
Tz(i,j)=Tz(i,j)-pt(k)*DetJac*k0*k0*dot(Et(i,:,k)',xxepsilonr(:,:,k)*Et(j,:,k)');
|
|
end
|
|
end
|
|
end
|
|
if domain(n)==5
|
|
for i=1:3
|
|
for k=1:IntegralOrder
|
|
Bz(i,1)= Bz(i,1)+pt(k)*DetJac*(k0*k0*dot(Et(i,:,k)',xxepsilonr(:,:,k)* Ei(:,k))-...
|
|
dot(curlEt(i,:,k)',xxinvXmur(:,:,k)*curlEi(:,k)));
|
|
end
|
|
end
|
|
end
|
|
|
|
|
|
for i=1:3
|
|
for j=1:3
|
|
A(EdgesOfElements((n-1)*3+i),EdgesOfElements((n-1)*3+j))=A(EdgesOfElements((n-1)*3+i),EdgesOfElements((n-1)*3+j))+Sz(i,j)+Tz(i,j);
|
|
end
|
|
b(EdgesOfElements((n-1)*3+i))=b(EdgesOfElements((n-1)*3+i))+Bz(i,1);
|
|
end
|
|
|
|
end
|
|
|
|
ii=0;
|
|
iii=0;
|
|
%边界积分
|
|
for n=1:NbrEdge*2
|
|
numEE=ElementNum(n,3);
|
|
if numEE==0
|
|
continue;
|
|
end
|
|
if ((edgeflag(numEE)==7)||(edgeflag(numEE)==8)||(edgeflag(numEE)==10)||(edgeflag(numEE)==11))
|
|
%坐标
|
|
ii=ii+1;
|
|
numEle=ElementNum(n,1);
|
|
numEdge=ElementNum(n,2);
|
|
x=zeros(5,1);y=zeros(5,1);l=zeros(3,1);
|
|
x(1)=node(Elements(1,numEle),1);y(1)=node(Elements(1,numEle),2);
|
|
x(2)=node(Elements(2,numEle),1);y(2)=node(Elements(2,numEle),2);
|
|
x(3)=node(Elements(3,numEle),1);y(3)=node(Elements(3,numEle),2);
|
|
x(4)=node(edge(1,numEE),1);y(4)=node(edge(1,numEE),2);
|
|
x(5)=node(edge(2,numEE),1);y(5)=node(edge(2,numEE),2);
|
|
l(1)=sqrt((x(1)-x(2))*(x(1)-x(2))+(y(1)-y(2))*(y(1)-y(2)));
|
|
l(2)=sqrt((x(1)-x(3))*(x(1)-x(3))+(y(1)-y(3))*(y(1)-y(3)));
|
|
l(3)=sqrt((x(3)-x(2))*(x(3)-x(2))+(y(3)-y(2))*(y(3)-y(2)));
|
|
|
|
avx=(x(1)+x(2)+x(3))/3;
|
|
avy=(y(1)+y(2)+y(3))/3;
|
|
rho2=avx*avx+avy*avy;
|
|
if rho2>1.5*1.5
|
|
iii=iii+1;
|
|
continue;
|
|
end
|
|
%雅可比矩阵
|
|
Jac=zeros(3,3);
|
|
Jac(1,1)=x(2)-x(1);Jac(1,2)=y(2)-y(1);
|
|
Jac(2,1)=x(3)-x(1);Jac(2,2)=y(3)-y(1);Jac(3,3)=1;
|
|
InvJac=inv(Jac);
|
|
%积分系数
|
|
phi1=(y(5) - y(4));phi2=(x(4) - x(5));
|
|
if x(5)~=x(4)
|
|
PhysicsFactor=(phi1)/(-phi2);
|
|
PhysicsFactor=abs(sqrt(1+PhysicsFactor*PhysicsFactor)*(-phi2)/2);
|
|
elseif x(5)==x(4)
|
|
PhysicsFactor=abs(phi1)/2;
|
|
end
|
|
%区域1外法向量
|
|
outerNormal=zeros(3,1);
|
|
outerNormal(1) = (phi1) / sqrt(phi1 * phi1 + phi2 * phi2);
|
|
outerNormal(2) = (phi2) / sqrt(phi1 * phi1 + phi2 * phi2);
|
|
%积分坐标修正
|
|
u=zeros(IntegralEdgeOrder,1);v=zeros(IntegralEdgeOrder,1);
|
|
if numEdge==1
|
|
for i=1:IntegralEdgeOrder
|
|
u(i)=(lx(i)+1)/2;
|
|
v(i)=0;
|
|
end
|
|
elseif numEdge==2
|
|
for i=1:IntegralEdgeOrder
|
|
u(i)=0;
|
|
v(i)=(lx(i)+1)/2;
|
|
end
|
|
elseif numEdge==3
|
|
for i=1:IntegralEdgeOrder
|
|
u(i)=(lx(i)+1)/2;
|
|
v(i)=(1-lx(i))/2;
|
|
end
|
|
end
|
|
%基函数
|
|
Et=zeros(3,3,IntegralEdgeOrder);curlEbt=zeros(3,IntegralEdgeOrder);temp=zeros(3,1);
|
|
for i=1:IntegralEdgeOrder
|
|
for j=1:3
|
|
[temp(1),temp(2)]=BF_Et(j,u(i),v(i));temp(3)=0;
|
|
Et(j,:,i)=InvJac*temp*l(j);
|
|
end
|
|
tempx=x(1)+Jac(1,1)*u(i)+Jac(2,1)*v(i);
|
|
tempy=y(1)+Jac(1,2)*u(i)+Jac(2,2)*v(i);
|
|
curlEbt(:,i)=CurlEinc(tempx,tempy);
|
|
end
|
|
|
|
%计算线积分矩阵
|
|
Bi=zeros(3,1);
|
|
for i=1:3
|
|
for k=1:IntegralEdgeOrder
|
|
Bi(i)=Bi(i)+lp(k)* PhysicsFactor*(dot(Et(i,:,k)',cross(-outerNormal,curlEbt(:,k))));
|
|
end
|
|
end
|
|
%矩阵组装
|
|
for i=1:3
|
|
b(EdgesOfElements((numEle-1)*3+i))=b(EdgesOfElements((numEle-1)*3+i))+Bi(i);
|
|
end
|
|
|
|
end
|
|
end
|
|
ii
|
|
iii
|
|
end |