163 lines
5.2 KiB
Matlab
163 lines
5.2 KiB
Matlab
function solver=assembly_inc(phy,mesh,solver)
|
|
|
|
k0=2*pi/phy.lda0;
|
|
%三角形面积分
|
|
[u,v,~,weight,nbrGP]=getGaussPoints(2);
|
|
w=1-u-v;
|
|
|
|
nbrOut=length(mesh.incIndex);
|
|
Ai=zeros(nbrOut*8*8,1);
|
|
Aj=zeros(nbrOut*8*8,1);
|
|
Av=zeros(nbrOut*8*8,1);
|
|
tempNN=0;
|
|
%loop of all tri
|
|
for n=1:nbrOut
|
|
domain=mesh.DomainOfTri(mesh.incIndex(n));
|
|
numTet=mesh.ConnOfTri(mesh.incIndex(n),1);
|
|
numFace=mesh.ConnOfTri(mesh.incIndex(n),2);
|
|
|
|
%vertex
|
|
x=mesh.Vertex(mesh.Tet(numTet,:),1);
|
|
y=mesh.Vertex(mesh.Tet(numTet,:),2);
|
|
z=mesh.Vertex(mesh.Tet(numTet,:),3);
|
|
|
|
MappingIndex=zeros(8,1);
|
|
if numFace==1 %123
|
|
x2=[1;0;0];
|
|
y2=[0;1;0];
|
|
z2=[0;0;1];
|
|
x3=[x(1);x(2);x(3)];
|
|
y3=[y(1);y(2);y(3)];
|
|
z3=[z(1);z(2);z(3)];
|
|
index=[1;3;7;...
|
|
2;4;8;...
|
|
13;14];
|
|
MappingIndex(1)=mesh.EdgeOfTet(numTet,1);
|
|
MappingIndex(2)=mesh.EdgeOfTet(numTet,2);
|
|
MappingIndex(3)=mesh.EdgeOfTet(numTet,4);
|
|
MappingIndex(4)=mesh.EdgeOfTet(numTet,1)+mesh.NbrEdge;
|
|
MappingIndex(5)=mesh.EdgeOfTet(numTet,2)+mesh.NbrEdge;
|
|
MappingIndex(6)=mesh.EdgeOfTet(numTet,4)+mesh.NbrEdge;
|
|
MappingIndex(7)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace);
|
|
MappingIndex(8)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace)+mesh.NbrFace;
|
|
elseif numFace==2 %124
|
|
x2=[1;0;0];
|
|
y2=[0;1;0];
|
|
z2=[0;0;0];
|
|
x3=[x(1);x(2);x(4)];
|
|
y3=[y(1);y(2);y(4)];
|
|
z3=[z(1);z(2);z(4)];
|
|
index=[1;5;9;...
|
|
2;6;10;...
|
|
15;16];
|
|
MappingIndex(1)=mesh.EdgeOfTet(numTet,1);
|
|
MappingIndex(2)=mesh.EdgeOfTet(numTet,3);
|
|
MappingIndex(3)=mesh.EdgeOfTet(numTet,5);
|
|
MappingIndex(4)=mesh.EdgeOfTet(numTet,1)+mesh.NbrEdge;
|
|
MappingIndex(5)=mesh.EdgeOfTet(numTet,3)+mesh.NbrEdge;
|
|
MappingIndex(6)=mesh.EdgeOfTet(numTet,5)+mesh.NbrEdge;
|
|
MappingIndex(7)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace);
|
|
MappingIndex(8)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace)+mesh.NbrFace;
|
|
elseif numFace==3 %134
|
|
x2=[1;0;0];
|
|
y2=[0;0;0];
|
|
z2=[0;1;0];
|
|
x3=[x(1);x(3);x(4)];
|
|
y3=[y(1);y(3);y(4)];
|
|
z3=[z(1);z(3);z(4)];
|
|
index=[3;5;11;...
|
|
4;6;12;...
|
|
17;18];
|
|
MappingIndex(1)=mesh.EdgeOfTet(numTet,2);
|
|
MappingIndex(2)=mesh.EdgeOfTet(numTet,3);
|
|
MappingIndex(3)=mesh.EdgeOfTet(numTet,6);
|
|
MappingIndex(4)=mesh.EdgeOfTet(numTet,2)+mesh.NbrEdge;
|
|
MappingIndex(5)=mesh.EdgeOfTet(numTet,3)+mesh.NbrEdge;
|
|
MappingIndex(6)=mesh.EdgeOfTet(numTet,6)+mesh.NbrEdge;
|
|
MappingIndex(7)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace);
|
|
MappingIndex(8)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace)+mesh.NbrFace;
|
|
elseif numFace==4 %234
|
|
x2=[0;0;0];
|
|
y2=[1;0;0];
|
|
z2=[0;1;0];
|
|
x3=[x(2);x(3);x(4)];
|
|
y3=[y(2);y(3);y(4)];
|
|
z3=[z(2);z(3);z(4)];
|
|
index=[7;9;11;...
|
|
8;10;12;...
|
|
19;20];
|
|
MappingIndex(1)=mesh.EdgeOfTet(numTet,4);
|
|
MappingIndex(2)=mesh.EdgeOfTet(numTet,5);
|
|
MappingIndex(3)=mesh.EdgeOfTet(numTet,6);
|
|
MappingIndex(4)=mesh.EdgeOfTet(numTet,4)+mesh.NbrEdge;
|
|
MappingIndex(5)=mesh.EdgeOfTet(numTet,5)+mesh.NbrEdge;
|
|
MappingIndex(6)=mesh.EdgeOfTet(numTet,6)+mesh.NbrEdge;
|
|
MappingIndex(7)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace);
|
|
MappingIndex(8)=mesh.NbrEdge*2+mesh.FaceOfTet(numTet,numFace)+mesh.NbrFace;
|
|
end
|
|
|
|
%jac
|
|
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);
|
|
|
|
%int coe
|
|
a=sqrt((x3(1)-x3(2))^2+(y3(1)-y3(2))^2+(z3(1)-z3(2))^2);
|
|
b=sqrt((x3(1)-x3(3))^2+(y3(1)-y3(3))^2+(z3(1)-z3(3))^2);
|
|
c=sqrt((x3(2)-x3(3))^2+(y3(2)-y3(3))^2+(z3(2)-z3(3))^2);
|
|
integCoe = 0.25 * sqrt((a + b + c) * (a + b - c) * (a - b + c) * (b + c - a));
|
|
|
|
%normal
|
|
normal=mesh.NormOfFace(domain,:)';
|
|
|
|
%bf
|
|
E=zeros(3,8,nbrGP); %3*基函数数目*GP节点数
|
|
for i=1:nbrGP
|
|
u2=x2(1)*u(i)+x2(2)*v(i)+x2(3)*w(i);
|
|
v2=y2(1)*u(i)+y2(2)*v(i)+y2(3)*w(i);
|
|
w2=z2(1)*u(i)+z2(2)*v(i)+z2(3)*w(i);
|
|
for j=1:8
|
|
E(:,j,i)=getBF(1,index(j),u2,v2,w2);
|
|
E(:,j,i)=Jac\E(:,j,i);
|
|
end
|
|
end
|
|
|
|
%mat
|
|
eps=phy.eps(mesh.DomainOfTet(numTet));
|
|
nn=sqrt(eps);
|
|
|
|
%submatrix
|
|
Ae=zeros(8,8);
|
|
for i=1:8
|
|
for j=1:8
|
|
for k=1:nbrGP
|
|
Ae(i,j)=Ae(i,j)+1i*k0*nn*integCoe*weight(k)*sum(E(:,i,k).*cross(normal,cross(E(:,j,k),normal)))*2;
|
|
end
|
|
end
|
|
end
|
|
Be=zeros(8,1);
|
|
for i=1:8
|
|
for k=1:nbrGP
|
|
Be(i)=Be(i)-1i*k0*nn*2*integCoe*weight(k)*sum(E(:,i,k).*cross(normal,cross(phy.Einc,normal)))*2;
|
|
end
|
|
end
|
|
|
|
%put in matrix
|
|
for i=1:8
|
|
for j=1:8
|
|
tempNN=tempNN+1;
|
|
Ai(tempNN)=MappingIndex(i);
|
|
Aj(tempNN)=MappingIndex(j);
|
|
Av(tempNN)=Ae(i,j);
|
|
end
|
|
solver.b(MappingIndex(i))=solver.b(MappingIndex(i))+Be(i);
|
|
end
|
|
|
|
end
|
|
|
|
solver.Ai=[solver.Ai;Ai];
|
|
solver.Aj=[solver.Aj;Aj];
|
|
solver.Av=[solver.Av;Av];
|
|
|
|
end |