You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
203 lines
4.9 KiB
203 lines
4.9 KiB
function smooth_compare_compliance()
|
|
nx=60;ny=1;nz=20;
|
|
load DH_afterSmooth25.mat DH
|
|
ke_all=localSM_all(nx,ny,nz,DH);
|
|
[U]=FEM(nx,ny,nz,ke_all);
|
|
comp_after = computeCompliance(nx,ny,nz,ke_all,U);
|
|
compliance_after=sum(sum(sum(comp_after)))
|
|
|
|
load DH_beforeSmooth25.mat DH
|
|
ke_all=localSM_all(nx,ny,nz,DH);
|
|
[U]=FEM(nx,ny,nz,ke_all);
|
|
comp_before = computeCompliance(nx,ny,nz,ke_all,U);
|
|
compliance_before=sum(sum(sum(comp_after)))
|
|
|
|
figure()
|
|
error=abs(comp_after-comp_before)./comp_before;
|
|
error=error.*(abs(comp_before)>0.2);
|
|
imagesc(squeeze(error(:,1,:)))
|
|
axis equal; axis tight; axis off;colorbar;
|
|
set(gca,'YDir','normal')
|
|
mean(error,'all')
|
|
end
|
|
|
|
function comp=computeCompliance(nx,ny,nz,ke_all,U)
|
|
comp =zeros(nz,ny,nx);
|
|
%situation 1,2
|
|
for elz = 1:nz
|
|
for ely = 1:ny
|
|
for elx = 1:nx
|
|
n2 = (nz+1)*(ny+1)*elx + (nz+1)*(ely-1) + elz;
|
|
n1 = n2-(nz+1)*(ny+1);
|
|
n3 = n1 + (nz+1);
|
|
n4 = n2 + (nz+1);
|
|
|
|
n5 = n1 + 1;
|
|
n6 = n2 + 1;
|
|
n7 = n3 + 1;
|
|
n8 = n4 + 1;
|
|
Ue = U([3*n1-2; 3*n1-1; 3*n1;
|
|
3*n2-2; 3*n2-1; 3*n2;
|
|
3*n3-2; 3*n3-1; 3*n3;
|
|
3*n4-2; 3*n4-1; 3*n4;
|
|
3*n5-2; 3*n5-1; 3*n5;
|
|
3*n6-2; 3*n6-1; 3*n6;
|
|
3*n7-2; 3*n7-1; 3*n7;
|
|
3*n8-2; 3*n8-1; 3*n8]);
|
|
|
|
comp(elz,ely,elx) = Ue'* ke_all(:,:,elx,ely,elz)*Ue;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
|
|
function [U]=FEM(nx,ny,nz,ke_all)
|
|
N = 3*(nx+1)*(ny+1)*(nz+1);
|
|
nele=nx*ny*nz;
|
|
si=zeros(1,nele*24^2);
|
|
sj=zeros(1,nele*24^2);
|
|
sv=zeros(1,nele*24^2);
|
|
|
|
for elx = 1:nx
|
|
for ely = 1:ny
|
|
for elz = 1:nz
|
|
n2 = (nz+1)*(ny+1)*elx + (nz+1)*(ely-1) + elz;
|
|
n1 = n2-(nz+1)*(ny+1);
|
|
n3 = n1 + (nz+1);
|
|
n4 = n2 + (nz+1);
|
|
n5 = n1 + 1;
|
|
n6 = n2 + 1;
|
|
n7 = n3 + 1;
|
|
n8 = n4 + 1;
|
|
edof = [3*n1-2; 3*n1-1; 3*n1;
|
|
3*n2-2; 3*n2-1; 3*n2;
|
|
3*n3-2; 3*n3-1; 3*n3;
|
|
3*n4-2; 3*n4-1; 3*n4;
|
|
3*n5-2; 3*n5-1; 3*n5;
|
|
3*n6-2; 3*n6-1; 3*n6;
|
|
3*n7-2; 3*n7-1; 3*n7;
|
|
3*n8-2; 3*n8-1; 3*n8];
|
|
ne=(elx-1)*ny*nz+(ely-1)*nz+elz;
|
|
si((ne-1)*24^2+1:ne*24^2)=kron(ones(1,24),edof.');
|
|
sj((ne-1)*24^2+1:ne*24^2)=kron(edof.',ones(1,24));
|
|
ke=ke_all(:,:,elx,ely,elz);
|
|
sv((ne-1)*24^2+1:ne*24^2)=ke(:);
|
|
% K(edof,edof) = K(edof,edof) + ke_all(:,:,elz,ely,elx);
|
|
end
|
|
end
|
|
end
|
|
K = sparse(si,sj,sv,N,N);
|
|
|
|
loadnodes=nx*(ny+1)*(nz+1)+1:nz+1:(nx+1)*(ny+1)*(nz+1)-nz;
|
|
loaddofs = loadnodes(:).*3;
|
|
F = zeros(N,1); F(loaddofs) = -1;
|
|
fixeddofs = [1:3*(ny+1)*(nz+1)];
|
|
alldofs = [1:N];
|
|
freedofs = setdiff(alldofs,fixeddofs);
|
|
U = zeros(N,1);
|
|
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);
|
|
end
|
|
|
|
|
|
function ke_all=localSM_all(nx,ny,nz,DH)
|
|
ke_all=zeros(24,24,nx,ny,nz);
|
|
xs=0;xd=1;
|
|
numgp_x=4;
|
|
gps_x=getgp(xs,xd,numgp_x);
|
|
weight_x=((xd-xs)/(numgp_x*2))^3;
|
|
for ix=1:numgp_x*2
|
|
x=gps_x(ix);
|
|
for iy=1:numgp_x*2
|
|
y=gps_x(iy);
|
|
for iz=1:numgp_x*2
|
|
z=gps_x(iz);
|
|
B=computeB(x,y,z);
|
|
for elx = 1:nx
|
|
for ely = 1:ny
|
|
for elz = 1:nz
|
|
Ce=DH{elx,ely,elz};
|
|
ke_all(:,:,elx,ely,elz)=ke_all(:,:,elx,ely,elz)+B.'*Ce*B;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
ke_all=ke_all.*weight_x;
|
|
end
|
|
|
|
function B=computeB(x,y,z)
|
|
B=zeros(6,24);
|
|
B(1,1)=-(y - 1)*(z - 1);
|
|
B(1,4)=(y - 1)*(z - 1);
|
|
B(1,7)=y*(z - 1);
|
|
B(1,10)=-y*(z - 1);
|
|
B(1,13)=z*(y - 1);
|
|
B(1,16)=-z*(y - 1);
|
|
B(1,19)=-y*z;
|
|
B(1,22)=y*z;
|
|
B(2,2)=-(x - 1)*(z - 1);
|
|
B(2,5)=x*(z - 1);
|
|
B(2,8)=(x - 1)*(z - 1);
|
|
B(2,11)=-x*(z - 1);
|
|
B(2,14)=z*(x - 1);
|
|
B(2,17)=-x*z;
|
|
B(2,20)=-z*(x - 1);
|
|
B(2,23)=x*z;
|
|
B(3,3)=-(x - 1)*(y - 1);
|
|
B(3,6)=x*(y - 1);
|
|
B(3,9)=y*(x - 1);
|
|
B(3,12)=-x*y;
|
|
B(3,15)=(x - 1)*(y - 1);
|
|
B(3,18)=-x*(y - 1);
|
|
B(3,21)=-y*(x - 1);
|
|
B(3,24)=x*y;
|
|
B(4,1)=-(x - 1)*(z - 1);
|
|
B(4,2)=-(y - 1)*(z - 1);
|
|
B(4,4)=x*(z - 1);
|
|
B(4,5)=(y - 1)*(z - 1);
|
|
B(4,7)=(x - 1)*(z - 1);
|
|
B(4,8)=y*(z - 1);
|
|
B(4,10)=-x*(z - 1);
|
|
B(4,11)=-y*(z - 1);
|
|
B(4,13)=z*(x - 1);
|
|
B(4,14)=z*(y - 1);
|
|
B(4,16)=-x*z;
|
|
B(4,17)=-z*(y - 1);
|
|
B(4,19)=-z*(x - 1);
|
|
B(4,20)=-y*z;
|
|
B(4,22)=x*z;
|
|
B(4,23)=y*z;
|
|
B(5,2)=-(x - 1)*(y - 1);
|
|
B(5,3)=-(x - 1)*(z - 1);
|
|
B(5,5)=x*(y - 1);
|
|
B(5,6)=x*(z - 1);
|
|
B(5,8)=y*(x - 1);
|
|
B(5,9)=(x - 1)*(z - 1);
|
|
B(5,11)=-x*y;
|
|
B(5,12)=-x*(z - 1);
|
|
B(5,14)=(x - 1)*(y - 1);
|
|
B(5,15)=z*(x - 1);
|
|
B(5,17)=-x*(y - 1);
|
|
B(5,18)=-x*z;
|
|
B(5,20)=-y*(x - 1);
|
|
B(5,21)=-z*(x - 1);
|
|
B(5,23)=x*y;
|
|
B(5,24)=x*z;
|
|
B(6,1)=-(x - 1)*(y - 1);
|
|
B(6,3)=-(y - 1)*(z - 1);
|
|
B(6,4)=x*(y - 1);
|
|
B(6,6)=(y - 1)*(z - 1);
|
|
B(6,7)=y*(x - 1);
|
|
B(6,9)=y*(z - 1);
|
|
B(6,10)=-x*y;
|
|
B(6,12)=-y*(z - 1);
|
|
B(6,13)=(x - 1)*(y - 1);
|
|
B(6,15)=z*(y - 1);
|
|
B(6,16)=-x*(y - 1);
|
|
B(6,18)=-z*(y - 1);
|
|
B(6,19)=-y*(x - 1);
|
|
B(6,21)=-y*z;
|
|
B(6,22)=x*y;
|
|
B(6,24)=y*z;
|
|
end
|