169 lines
5.3 KiB
Matlab
169 lines
5.3 KiB
Matlab
function [A,b,EdgesOfElements]=FemAssemble3(node,Elements,domain,lam0,epsilonr,mur)
|
||
%圆形PML矩阵组装
|
||
|
||
%物理参数
|
||
k0=2*pi/lam0;
|
||
E0=[0,1];
|
||
%高斯积分点
|
||
[xt,yt,pt,IntegralOrder]=GetGuassPoints(4);
|
||
% IntegralOrder=4;
|
||
% xt= [ 0.333333333333333,0.6,0.2,0.2 ];
|
||
% yt= [ 0.3333333333333333,0.2,0.6,0.2 ];
|
||
% pt= [ -0.28125,.260416666666,.260416666666,.260416666666 ];
|
||
%网格处理
|
||
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));
|
||
%矩阵初始化
|
||
Dof=NbrEdges;
|
||
A=zeros(Dof,Dof);
|
||
b=zeros(Dof,1);
|
||
%狄利克雷边界
|
||
NbrDboudary=0;
|
||
Dboudary=zeros(NbrEdges,2);
|
||
|
||
%循环网格 组装矩阵
|
||
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
|
||
|
||
|
||
xxepsilonr=zeros(3,3,IntegralOrder);
|
||
xxinvXmur=zeros(3,3,IntegralOrder);
|
||
if domain(n)==4
|
||
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=100;
|
||
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);
|
||
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
|
||
|
||
|
||
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
|
||
end
|
||
|
||
%边界处理
|
||
%NbrDboudary 边界数
|
||
%Dboudary n*2 1:网格 2：第几条边
|
||
if (x(1)==0)&&(x(2)==0&&(abs(y(1))<1.01)&&(abs(y(2))<1.01))
|
||
NbrDboudary=NbrDboudary+1;
|
||
Dboudary(NbrDboudary,1)=n;
|
||
Dboudary(NbrDboudary,2)=1;
|
||
elseif (x(1)==0)&&(x(3)==0&&(abs(y(1))<1.01)&&(abs(y(3))<1.01))
|
||
NbrDboudary=NbrDboudary+1;
|
||
Dboudary(NbrDboudary,1)=n;
|
||
Dboudary(NbrDboudary,2)=2;
|
||
elseif (x(2)==0)&&(x(3)==0&&(abs(y(2))<1.01)&&(abs(y(3))<1.01))
|
||
NbrDboudary=NbrDboudary+1;
|
||
Dboudary(NbrDboudary,1)=n;
|
||
Dboudary(NbrDboudary,2)=3;
|
||
end
|
||
|
||
end
|
||
|
||
|
||
%强加狄利克雷边界条件
|
||
for n=1:NbrDboudary
|
||
numEle=Dboudary(n,1);
|
||
outNorm=[-1;0];
|
||
if (Dboudary(n,2)==1)
|
||
u=0.5;v=0;
|
||
elseif (Dboudary(n,2)==2)
|
||
u=0;v=0.5;
|
||
elseif (Dboudary(n,2)==3)
|
||
u=0.5;v=0.5;
|
||
end
|
||
%节点坐标
|
||
x=zeros(3,1);y=zeros(3,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);
|
||
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);
|
||
[temp(1),temp(2)]=BF_Et(Dboudary(n,2),u,v);temp(3)=0;
|
||
tempEt=InvJac*temp*l(Dboudary(n,2));
|
||
temp2=outNorm(1)*tempEt(2)-outNorm(2)*tempEt(1);
|
||
|
||
numEdge=EdgesOfElements((Dboudary(n,1)-1)*3+Dboudary(n,2));
|
||
temp1=outNorm(1)*E0(2)-outNorm(2)*E0(1);
|
||
pp=temp1/temp2;
|
||
b(numEdge)=pp;
|
||
A(numEdge,:)=0;
|
||
b(:)=b(:)-pp*A(:,numEdge);
|
||
A(:,numEdge)=0;
|
||
A(numEdge,numEdge)=1;
|
||
end
|
||
|
||
end |