function [Kgl, M] = K_assembly (in_data, type_analysis, part_N_r, part_N_c)
[dofPnd,EL_TYPE] = type_of_elem (in_data);
dof = size(in_data.ND)*dofPnd;
nn_z = floor((1+part_N_r(2)-part_N_r(1))*(1+part_N_c(2)-part_N_c(1))*0.05);
Kgl = spalloc(1+part_N_r(2)-part_N_r(1), 1+part_N_c(2)-part_N_c(1),nn_z);
if isfield(in_data,'T')
M = spalloc(1+part_N_r(2)-part_N_r(1), 1,nn_z);
else
M = spalloc(1+part_N_r(2)-part_N_r(1), 1+part_N_c(2)-part_N_c(1),nn_z);
end
indj =1;
sxEL = size(in_data.EL,1);
for i=1:sxEL
if i>(indj*250)
disp([' K assembly: ' num2str(indj*250) ' elements ']);
indj=indj+1;
end;
if i==sxEL
disp([' K assembly: ' num2str(i) ' elements - complete']);
end;
if in_data.EL(i,2)==44
dofPnd = 1;
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
a=in_data.T.ND;
in_data.T.Tmatrix = zeros(size(in_data.EL,1),3);
in_data.T.hmatrix = zeros(size(in_data.EL,1),3);
for i=1:size(in_data.EL,1)
b1 = in_data.EL(i,3:4);
b2 = in_data.EL(i,4:5);
b3 = [in_data.EL(i,5) in_data.EL(i,3)];
ind1=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(1) ones(1,length(a)-1)'*b1(2)]')));
ind1_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(2) ones(1,length(a)-1)'*b1(1)]')));
ind2=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(1) ones(1,length(a)-1)'*b2(2)]')));
ind2_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(2) ones(1,length(a)-1)'*b2(1)]')));
ind3=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(1) ones(1,length(a)-1)'*b3(2)]')));
ind3_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(2) ones(1,length(a)-1)'*b3(1)]')));
if [ind1 ind1_]
in_data.T.Tmatrix(i,1) = in_data.T.T([ind1 ind1_]);
in_data.T.hmatrix(i,1)= in_data.T.h([ind1 ind1_]);
end
if [ind2 ind2_]
in_data.T.Tmatrix(i,2) = in_data.T.T([ind2 ind2_]);
in_data.T.hmatrix(i,2)= in_data.T.h([ind2 ind2_]);
end
if [ind3 ind3_]
in_data.T.Tmatrix(i,3) = in_data.T.T([ind3 ind3_]);
in_data.T.hmatrix(i,3)= in_data.T.h([ind3 ind3_]);
end
end
a=in_data.flux.ND;
in_data.flux.qmatrix = zeros(size(in_data.EL,1),3);
for i=1:size(in_data.EL,1)
b1 = in_data.EL(i,3:4);
b2 = in_data.EL(i,4:5);
b3 = [in_data.EL(i,5) in_data.EL(i,3)];
ind1=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(1) ones(1,length(a)-1)'*b1(2)]')));
ind1_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(2) ones(1,length(a)-1)'*b1(1)]')));
ind2=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(1) ones(1,length(a)-1)'*b2(2)]')));
ind2_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(2) ones(1,length(a)-1)'*b2(1)]')));
ind3=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(1) ones(1,length(a)-1)'*b3(2)]')));
ind3_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(2) ones(1,length(a)-1)'*b3(1)]')));
if [ind1 ind1_]
in_data.flux.qmatrix(i,1) = in_data.flux.q([ind1 ind1_]);
end
if [ind2 ind2_]
in_data.flux.qmatrix(i,2) = in_data.flux.q([ind2 ind2_]);
end
if [ind3 ind3_]
in_data.flux.qmatrix(i,3) = in_data.flux.q([ind3 ind3_]);
end
end
kx = in_data.EL(i,6);
ky = in_data.EL(i,7);
Q = in_data.EL(i,8);
h = in_data.T.hmatrix(i,:);
Tinf = in_data.T.Tmatrix(i,:);
q = in_data.flux.qmatrix(i,:);
[Klc,rlc] = D2_CST_temp (in_data.ND(node1,2),in_data.ND(node1,3), ...
in_data.ND(node2,2),in_data.ND(node2,3), ...
in_data.ND(node3,2),in_data.ND(node3,3), ...
kx, ky, h, Q, Tinf, q);
t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd;
bg = [(t1-(dofPnd-1)):t1]; md = [(t2-(dofPnd-1)):t2];
en = [(t3-(dofPnd-1)):t3];
Kgl([bg md en],[bg md en]) = Kgl([bg md en],[bg md en]) + Klc;
M([bg md en]) = M([bg md en]) + rlc;
end
if in_data.EL(i,2)==0 | in_data.EL(i,2)==1 | in_data.EL(i,2)==2 | in_data.EL(i,2)==31
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,5); end;
if in_data.mater.A~=0 A_1=in_data.mater.A; else A_1 = in_data.EL(i,6); end;
if in_data.mater.I~=0 I_1=in_data.mater.I; else I_1 = in_data.EL(i,7); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_= in_data.EL(i,8); end;
if in_data.EL(i,2)==0
[M_loc,Klc] = FF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),...
in_data.ND(node2,2), in_data.ND(node2,3),E1,A_1,I_1,rho_);
end;
if in_data.EL(i,2)==1
[M_loc,Klc] = FP_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),...
in_data.ND(node2,2), in_data.ND(node2,3),E1,A_1,I_1,rho_);
end;
if in_data.EL(i,2)==2
[M_loc,Klc] = PF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),...
in_data.ND(node2,2), in_data.ND(node2,3),E1,A_1,I_1,rho_);
end;
if in_data.EL(i,2)==31
[M_loc,Klc] = D3_TRUSS (in_data.ND(node1,2),in_data.ND(node1,3),...
in_data.ND(node1,4), in_data.ND(node2,2), in_data.ND(node2,3),...
in_data.ND(node2,4), E1, A_1, rho_);
end;
t1 = node1*dofPnd; t2 = node2*dofPnd;
bg = [t1-(dofPnd-1):t1]; en = [t2-(dofPnd-1):t2];
Kgl([bg en],[bg en]) = Kgl([bg en],[bg en])+Klc;
M([bg en],[bg en]) = M([bg en],[bg en])+M_loc;
end;
if in_data.EL(i,2)==4
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,6); end;
if in_data.mater.h~=0 h_1=in_data.mater.h; else h_1 = in_data.EL(i,7); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,8); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_= in_data.EL(i,9); end;
[Bsys,Esys,Klc,Msys] = D2_CST (in_data.ND(node1,2),in_data.ND(node1,3), ...
in_data.ND(node2,2),in_data.ND(node2,3), in_data.ND(node3,2),...
in_data.ND(node3,3), E1,h_1,miu_1);
Msys = Msys.*rho_;
t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd;
bg = [(t1-(dofPnd-1)):t1]; md = [(t2-(dofPnd-1)):t2];
en = [(t3-(dofPnd-1)):t3];
Kgl([bg md en],[bg md en]) = Kgl([bg md en],[bg md en]) + Klc;
M([bg md en],[bg md en]) = M([bg md en],[bg md en]) + Msys;
end;
if in_data.EL(i,2)==3
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,5); end;
if in_data.mater.G~=0 G1=in_data.mater.G; else G1 = in_data.EL(i,6); end;
if in_data.mater.b~=0 b1=in_data.mater.b; else b1 = in_data.EL(i,7); end;
if in_data.mater.ho~=0 ho1=in_data.mater.ho; else ho1 = in_data.EL(i,8); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_= in_data.EL(i,9); end;
[M_loc,Klc] = D3_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node1,4), ...
in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node2,4),E1,G1,b1,ho1,rho_);
t1 = node1*dofPnd; t2 = node2*dofPnd;
bg = [(t1-(dofPnd-1)):t1]; en = [(t2-(dofPnd-1)):t2];
Kgl([bg en],[bg en]) = Kgl([bg en],[bg en]) + Klc;
M([bg en],[bg en]) = M([bg en],[bg en]) + M_loc;
end;
if in_data.EL(i,2)==6
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
node5 = find(in_data.ND(:,1)==in_data.EL(i,7));
node6 = find(in_data.ND(:,1)==in_data.EL(i,8));
node7 = find(in_data.ND(:,1)==in_data.EL(i,9));
node8 = find(in_data.ND(:,1)==in_data.EL(i,10));
x1 = in_data.ND(node1,2); y1 = in_data.ND(node1,3); z1 = in_data.ND(node1,4);
x2 = in_data.ND(node2,2); y2 = in_data.ND(node2,3); z2 = in_data.ND(node2,4);
x3 = in_data.ND(node3,2); y3 = in_data.ND(node3,3); z3 = in_data.ND(node3,4);
x4 = in_data.ND(node4,2); y4 = in_data.ND(node4,3); z4 = in_data.ND(node4,4);
x5 = in_data.ND(node5,2); y5 = in_data.ND(node5,3); z5 = in_data.ND(node5,4);
x6 = in_data.ND(node6,2); y6 = in_data.ND(node6,3); z6 = in_data.ND(node6,4);
x7 = in_data.ND(node7,2); y7 = in_data.ND(node7,3); z7 = in_data.ND(node7,4);
x8 = in_data.ND(node8,2); y8 = in_data.ND(node8,3); z8 = in_data.ND(node8,4);
if in_data.mater.E~=0 Em=in_data.mater.E; else Em=in_data.EL(i,11); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,12); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_=in_data.EL(i,13); end;
[Klc,Bsys,Esys,V] = D3_LB (x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,...
x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,Em,miu_1);
sf = ones(24,1);
fM = V*rho_/8.*sf';
t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd; t4 = node4*dofPnd;
t5 = node5*dofPnd; t6 = node6*dofPnd; t7 = node7*dofPnd; t8 = node8*dofPnd;
bg = [(t1-(dofPnd-1)):t1];
md1= [(t2-(dofPnd-1)):t2];
md2= [(t3-(dofPnd-1)):t3];
md3= [(t4-(dofPnd-1)):t4];
md4= [(t5-(dofPnd-1)):t5];
md5= [(t6-(dofPnd-1)):t6];
md6= [(t7-(dofPnd-1)):t7];
en = [(t8-(dofPnd-1)):t8];
Kgl([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) = ...
Kgl([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) + Klc;
M([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) = ...
M([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) + diag(fM);
end;
if in_data.EL(i,2)==10
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
if in_data.mater.E~=0 Em=in_data.mater.E; else Em=in_data.EL(i,7); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,8); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_=in_data.EL(i,9); end;
[Klc,Mlc] = D3_TETRAH (in_data.ND(node1,2), in_data.ND(node1,3), in_data.ND(node1,4),...
in_data.ND(node2,2), in_data.ND(node2,3), in_data.ND(node2,4),...
in_data.ND(node3,2), in_data.ND(node3,3), in_data.ND(node3,4),...
in_data.ND(node4,2), in_data.ND(node4,3), in_data.ND(node4,4), Em, miu_1);
Mlc = Mlc.*rho_;
t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd; t4 = node4*dofPnd;
bg = [(t1-(dofPnd-1)):t1]; md1= [(t2-(dofPnd-1)):t2];
md2= [(t3-(dofPnd-1)):t3]; md3= [(t4-(dofPnd-1)):t4];
Kgl([bg md1 md2 md3],[bg md1 md2 md3]) = ...
Kgl([bg md1 md2 md3],[bg md1 md2 md3]) + Klc;
M([bg md1 md2 md3],[bg md1 md2 md3]) = ...
M([bg md1 md2 md3],[bg md1 md2 md3]) + Mlc;
end;
if in_data.EL(i,2)==5
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
if in_data.mater.E~=0 Em=in_data.mater.E; else Em=in_data.EL(i,7); end;
if in_data.mater.h~=0
if length(in_data.mater.h)==1 hT=[1 1 1 1].*in_data.mater.h; end
if length(in_data.mater.h)==4 hT=in_data.mater.h; end
else hT = in_data.EL(i,8:11); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,12); end;
if in_data.mater.rho~=0 rho_=in_data.mater.rho; else rho_=in_data.EL(i,13); end;
if in_data.EL(i,2)==5
[Bsys,Esys,Klc,Msys] = D2_CSQ (in_data.ND(node1,2),in_data.ND(node1,3), ...
in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node3,2),in_data.ND(node3,3),...
in_data.ND(node4,2),in_data.ND(node4,3), Em, hT, miu_1);
end
if in_data.EL(i,2)==51
ex = in_data.PML.ex;
ey = in_data.PML.ey;
[Bsys,Esys,Klc,Msys] = D2_CSQ_PML (in_data.ND(node1,2),in_data.ND(node1,3), ...
in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node3,2),in_data.ND(node3,3),...
in_data.ND(node4,2),in_data.ND(node4,3), Em, hT, miu_1, ex, ey);
end
M_loc = spalloc(8,8,1);
for z=1:8 M_loc(z,z)=1; end;
V_loc = .5*(( in_data.ND(node3,2)-in_data.ND(node1,2) )*...
( in_data.ND(node4,3)-in_data.ND(node2,3) )-...
( in_data.ND(node4,2)-in_data.ND(node2,2) )*...
( in_data.ND(node3,3)-in_data.ND(node1,3) ))* mean(hT);
Msys = M_loc.*rho_*V_loc/4;
if in_data.EL(i,2)==51
Msys = Msys.*diag([ex ey ex ey ex ey ex ey]);
end
t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd; t4 = node4*dofPnd;
bg = [(t1-(dofPnd-1)):t1]; md1= [(t2-(dofPnd-1)):t2];
md2= [(t3-(dofPnd-1)):t3]; en = [(t4-(dofPnd-1)):t4];
Kgl([bg md1 md2 en],[bg md1 md2 en]) = ...
Kgl([bg md1 md2 en],[bg md1 md2 en]) + Klc;
M([bg md1 md2 en],[bg md1 md2 en]) = ...
M([bg md1 md2 en],[bg md1 md2 en]) + Msys;
end;
if in_data.EL(i,2)==9
node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
if in_data.mater.E~=0 E1=in_data.mater.E; else E1=in_data.EL(i,6); end;
if in_data.mater.h~=0 h_1=in_data.mater.h; else h_1 = in_data.EL(i,7); end;
if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else miu_1 = in_data.EL(i,8); end;
x1 = in_data.ND(node1,2); y1 = in_data.ND(node1,3);
x2 = in_data.ND(node2,2); y2 = in_data.ND(node2,3);
x3 = in_data.ND(node3,2); y3 = in_data.ND(node3,3);
[Klc,ML] = D2_BCIZ (x1,y1,x2,y2,x3,y3,E1,h_1,miu_1);
t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd;
A = (0.5 * det([1 1 1; x1 x2 x3; y1 y2 y3]) );
ML = ML.*in_data.mater.rhoX;
bg = [(t1-2):t1]; md1= [(t2-2):t2]; md2= [(t3-2):t3];
Kgl([bg md1 md2],[bg md1 md2]) = Kgl([bg md1 md2],[bg md1 md2]) + Klc;
M([bg md1 md2],[bg md1 md2]) = M([bg md1 md2],[bg md1 md2]) + ML;
end;
end;
sp_=sprintf('\n');
disp([sp_ ' nnz in K: ' num2str(nnz(Kgl))]);
