Wxwei
/
MYPINN
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
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.
5 Commits
1 Branch
0 Tags
72 MiB
 
 
 
Branch: master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
MYPINN/numInvHom_code/inverseHomog3d.m

341 lines
14 KiB
Raw Permalink Blame History

% Q: the elastic tensor
% dQ: the derivatives of elastic tensor for microstructures
function [Q,dQ,xPhys] = inverseHomog3d()
E = 0.1; nu = 0.15;
% E = 0.1; nu = 0.35;
% E = 0.1; nu = 0.05;
D_target = E/(1+nu)/(1-2*nu)* [ 1-nu nu nu 0 0 0;
nu 1-nu nu 0 0 0;
nu nu 1-nu 0 0 0;
0 0 0 (1-2*nu)/2 0 0;
0 0 0 0 (1-2*nu)/2 0;
0 0 0 0 0 (1-2*nu)/2;];
nelx = 20;
nely = 20;
nelz = 20;
nele = nelx*nely*nelz;
volfrac = 0.5;
rmin = 2;
E0 = 1;
Emin = 1e-9;
nu = 0.3;
penal = 5;
% the initial definitions of the PUC
D0 = E0/(1+nu)/(1-2*nu)* [ 1-nu nu nu 0 0 0;
nu 1-nu nu 0 0 0;
nu nu 1-nu 0 0 0;
0 0 0 (1-2*nu)/2 0 0;
0 0 0 0 (1-2*nu)/2 0;
0 0 0 0 0 (1-2*nu)/2;];
% dx = lx/nelx; dy = ly/nely; dz = lz/nelz;
% KE = elementMatVec3D(1, 1, 1, D0);
KE = elementMatVec3D(1.0/2, 1.0/2, 1.0/2, D0);
% [KE] = lk_H8(E0,nu);
Num_node = (1+nely)*(1+nelx)*(1+nelz);
nodenrs = reshape(1:Num_node,1+nely,1+nelx,1+nelz);
edofVec = reshape(3*nodenrs(1:end-1,1:end-1,1:end-1)+1,nelx*nely*nelz,1);
edofMat = repmat(edofVec,1,24)+repmat([0 1 2 3*nely+[3 4 5 0 1 2] -3 -2 -1 3*(nelx+1)*(nely+1)+[0 1 2 3*nely+[3 4 5 0 1 2] -3 -2 -1]], nele, 1);
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
[H,Hs] = filter(nelx,nely,nelz,rmin);
% periodic boundary conditions
n1 = [nodenrs(end, [1 end], 1) nodenrs(1, [end 1], 1) nodenrs(end, [1 end], end) nodenrs(1, [end 1], end)];
d1 = reshape([3*n1-2; 3*n1-1; 3*n1],3*numel(n1),1);
n3 = [reshape(squeeze(nodenrs(end,1,2:end-1)),1,numel(squeeze(nodenrs(end,1,2:end-1))))...
reshape(squeeze(nodenrs(1, 1, 2:end-1)),1,numel(squeeze(nodenrs(1, 1, 2:end-1))))...
reshape(squeeze(nodenrs(end,2:end-1,1)),1,numel(squeeze(nodenrs(end,2:end-1,1))))...
reshape(squeeze(nodenrs(1, 2:end-1, 1)),1,numel(squeeze(nodenrs(1, 2:end-1, 1))))...
reshape(squeeze(nodenrs(2:end-1, 1, 1)),1,numel(squeeze(nodenrs(2:end-1, 1, 1))))...
reshape(squeeze(nodenrs(2:end-1,1,end)),1,numel(squeeze(nodenrs(2:end-1,1,end))))...
reshape(squeeze(nodenrs(2:end-1, 2:end-1, 1)),1,numel(squeeze(nodenrs(2:end-1, 2:end-1, 1))))...
reshape(squeeze(nodenrs(2:end-1, 1, 2:end-1)),1,numel(squeeze(nodenrs(2:end-1, 1, 2:end-1))))...
reshape(squeeze(nodenrs(end,2:end-1,2:end-1)),1,numel(squeeze(nodenrs(end,2:end-1,2:end-1))))];
d3 = reshape([3*n3-2; 3*n3-1; 3*n3],3*numel(n3),1);
n4 = [reshape(squeeze(nodenrs(1, end, 2:end-1)),1,numel(squeeze(nodenrs(1, end, 2:end-1))))...
reshape(squeeze(nodenrs(end,end,2:end-1)),1,numel(squeeze(nodenrs(end,end,2:end-1))))...
reshape(squeeze(nodenrs(1, 2:end-1, end)),1,numel(squeeze(nodenrs(1, 2:end-1, end))))...
reshape(squeeze(nodenrs(end,2:end-1,end)),1,numel(squeeze(nodenrs(end,2:end-1,end))))...
reshape(squeeze(nodenrs(2:end-1,end,end)),1,numel(squeeze(nodenrs(2:end-1,end,end))))...
reshape(squeeze(nodenrs(2:end-1, end, 1)),1,numel(squeeze(nodenrs(2:end-1, end, 1))))...
reshape(squeeze(nodenrs(2:end-1,2:end-1,end)),1,numel(squeeze(nodenrs(2:end-1,2:end-1,end))))...
reshape(squeeze(nodenrs(2:end-1,end,2:end-1)),1,numel(squeeze(nodenrs(2:end-1,end,2:end-1))))...
reshape(squeeze(nodenrs(1, 2:end-1, 2:end-1)),1,numel(squeeze(nodenrs(1, 2:end-1, 2:end-1))))];
d4 = reshape([3*n4-2; 3*n4-1; 3*n4],3*numel(n4),1);
n2 = setdiff(nodenrs(:),[n1(:);n3(:);n4(:)]);
d2 = reshape([3*n2-2; 3*n2- 1; 3*n2],3*numel(n2),1);
% the imposing of six linearly independent unit test strains
e = eye(6); ufixed = zeros(24,6);
vert_cor = [0, nelx, nelx, 0, 0, nelx, nelx, 0;
0, 0, nely, nely, 0, 0, nely, nely;
0, 0, 0, 0, nelz, nelz, nelz, nelz];
for i = 1:6
epsilon = [ e(i,1), e(i,4)/2, e(i,6)/2;
e(i,4)/2, e(i,2), e(i,5)/2;
e(i,6)/2, e(i,5)/2, e(i,3)];
ufixed(:,i) = reshape(epsilon*vert_cor,24,1);
end
% 3D boundary constraint equations
wfixed = [repmat(ufixed(7:9,:),numel(squeeze(nodenrs(end,1,2:end-1))),1);
repmat(ufixed( 4:6,:)-ufixed(10:12,:),numel(squeeze(nodenrs(1, 1, 2:end-1))),1);
repmat(ufixed(22:24,:),numel(squeeze(nodenrs(end,2:end-1,1))),1);
repmat(ufixed(13:15,:)-ufixed(10:12,:),numel(squeeze(nodenrs(1, 2:end-1, 1))),1);
repmat(ufixed(16:18,:),numel(squeeze(nodenrs(2:end-1, 1, 1))),1);
repmat(ufixed( 4:6,:)-ufixed(13:15,:),numel(squeeze(nodenrs(2:end-1,1,end))),1);
repmat(ufixed(13:15,:),numel(squeeze(nodenrs(2:end-1, 2:end-1,1))),1);
repmat(ufixed(4:6,:),numel(squeeze(nodenrs(2:end-1, 1, 2:end- 1))),1);
repmat(ufixed(10:12,:),numel(squeeze(nodenrs(end,2:end-1,2:end-1))),1)];
% intial
qe = cell(6,6);
Q = zeros(6,6);
dQ = cell(6,6);
x = ones(nely,nelx,nelz).*volfrac;
x(nely/2:nely/2+1,nelx/2:nelx/2+1,nelz /2:nelz/2+1) = Emin;
% x = max(Emin,rand(nely,nelx,nelz));
xPhys = x;
% beta = 1;
% xTilde = x;
% xPhys = 1-exp(-beta*xTilde)+xTilde*exp(-beta);
loopbeta = 0; loop = 0; maxloop = 100;
change = 1;
% figure;
OBJ = zeros(maxloop,1);
VOL = zeros(maxloop,1);
BWR = zeros(maxloop,1);
while loop < maxloop && change > 1e-4
loop = loop+1; loopbeta = loopbeta+1;
% FE-ANALYSIS
sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(1-Emin)),24*24*nele,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
Kr = [K(d2,d2), K(d2,d3)+K(d2,d4);
K(d3,d2)+K(d4,d2),K(d3,d3)+K(d4,d3)+K(d3,d4)+K(d4,d4)];
U(d1,:) = ufixed;
U([d2;d3],:) = Kr\(-[K(d2,d1);K(d3,d1)+K(d4,d1)]*ufixed-[K(d2,d4);K(d3,d4)+K(d4,d4)]*wfixed);
% for i_load=1:6
% U([d2;d3],i_load) = pcg(Kr, (-[K(d2,d1);K(d3,d1)+K(d4,d1)]*ufixed(:,i_load)-[K(d2,d4);K(d3,d4)+K(d4,d4)]*wfixed(:,i_load)),1e-3,200);
% end
U(d4,:) = U(d3,:)+wfixed;
% homog
cellVolume = nele;
for i = 1:6
for j = 1:6
U1 = U(:,i);
U2 = U(:,j);
qe{i,j} = reshape(sum((U1(edofMat)*KE).*U2(edofMat),2),nely,nelx,nelz);
Q(i,j) = 1/cellVolume*sum(sum(sum((Emin+xPhys.^penal*(1-Emin)).*qe{i,j})));
dQ{i,j} = 1/cellVolume*(penal*(1-Emin)*xPhys.^(penal-1).*qe{i,j});
end
end
% max shear modulus
% c = -(Q(4,4)+Q(5,5)+Q(6,6));
% dc = -(dQ{4,4}+dQ{5,5}+dQ{6,6});
% dQ_ele = zeros(nely,nelx,nelz,6,6);
% for i=1:6
% for j=1:6
% dQ_ele(:,:,:,i,j) = dQ{i,j};
% end
% end
% c = sum((Q-D_target).^2,'all');
% dc = zeros(nely,nelx,nelz);
% for ix = 1:nelx
% for iy = 1:nely
% for iz = 1:nelz
% dc(iy,ix,iz) = sum((Q-D_target).*squeeze(dQ_ele(iy,ix,iz,:,:)).*2,'all');
% end
% end
% end
% [c,dc] = matchCH2(Q,D_target,'frob2', dQ);
c = sum((Q-D_target).^2,'all');
dQ_ele = zeros(nely,nelx,nelz,6,6);
for i=1:6
for j=1:6
dQ_ele(:,:,:,i,j) = dQ{i,j};
end
end
dc = zeros(nely,nelx,nelz);
for ix = 1:nelx
for iy = 1:nely
for iz = 1:nelz
dc(iy,ix,iz) = sum((Q-D_target).'.*squeeze(dQ_ele(iy,ix,iz,:,:)).*2,'all');
end
end
end
disp([min(dc,[],'all'), max(dc,[],'all')])
dv = ones(nely,nelx,nelz);
% FILTERING AND MODIFICATION OF SENSITIVITIES
dc(:) = H*(dc(:)./Hs);
dv(:) = H*(dv(:)./Hs);
% OPTIMALITY CRITERIA UPDATE
l1 = 0; l2 = 1e9; move = 0.2;
while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1);
% xnew = max(Emin,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
xnew = max(Emin,max(x-move,min(1,min(x+move,x.*(-dc./dv/lmid)))));
xPhys(:) = (H*xnew(:))./Hs;
if sum(xPhys(:)) > volfrac*nele, l1 = lmid; else l2 = lmid; end
end
change = max(abs(xnew(:)-x(:)));
x = xnew;
% PRINT RESULTS
fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c,mean(xPhys(:)),change);
disp(['black-white ratio : ',num2str(max(xPhys,[],'all')-min(xPhys,[],'all'))])
% PLOT DENSITIES
% clf; display_3D(xPhys);
OBJ(loop) = c;
VOL(loop) = mean(xPhys(:));
BWR(loop) = max(xPhys,[],'all')-min(xPhys,[],'all');
plot_convergence_all(loop,maxloop,OBJ,VOL,BWR);
end
figure();clf; display_3D(xPhys);
end
% === DISPLAY 3D TOPOLOGY (ISO-VIEW) ===
function display_3D(rho)
[nely,nelx,nelz] = size(rho);
hx = 1; hy = 1; hz = 1; % User-defined unit element size
face = [1 2 3 4; 2 6 7 3; 4 3 7 8; 1 5 8 4; 1 2 6 5; 5 6 7 8];
set(gcf,'Name','ISO display','NumberTitle','off');
for k = 1:nelz
z = (k-1)*hz;
for i = 1:nelx
x = (i-1)*hx;
for j = 1:nely
y = nely*hy - (j-1)*hy;
if (rho(j,i,k) > 0.5) % User-defined display density threshold
vert = [x y z; x y-hx z; x+hx y-hx z; x+hx y z; x y z+hx;x y-hx z+hx; x+hx y-hx z+hx;x+hx y z+hx];
vert(:,[2 3]) = vert(:,[3 2]); vert(:,2,:) = -vert(:,2,:);
patch('Faces',face,'Vertices',vert,'FaceColor',[0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k))]);
hold on;
end
end
end
end
axis equal; axis tight; axis off; box on; view([30,30]); pause(1e-6);
end
% PREPARE FILTER
function [H,Hs] = filter(nelx,nely,nelz,rmin)
nele = nelx*nely*nelz;
iH = ones(nele*(2*(ceil(rmin)-1)+1)^2,1);
jH = ones(size(iH));
sH = zeros(size(iH));
k = 0;
for k1 = 1:nelz
for i1 = 1:nelx
for j1 = 1:nely
e1 = (k1-1)*nelx*nely + (i1-1)*nely+j1;
for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz)
for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
e2 = (k2-1)*nelx*nely + (i2-1)*nely+j2;
k = k+1;
iH(k) = e1;
jH(k) = e2;
sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2+(k1-k2)^2));
end
end
end
end
end
end
H = sparse(iH,jH,sH);
Hs = sum(H,2);
end
function Ke = elementMatVec3D(a, b, c, DH)
GN_x=[-1/sqrt(3),1/sqrt(3)];
GN_y=GN_x;
GN_z=GN_x;
GaussWeigh=[1,1];
Ke = zeros(24,24); L = zeros(6,9);
L(1,1) = 1; L(2,5) = 1; L(3,9) = 1;
L(4,2) = 1; L(4,4) = 1; L(5,6) = 1;
L(5,8) = 1; L(6,3) = 1; L(6,7) = 1;
for i=1:length(GN_x)
for j=1:length(GN_y)
for k=1:length(GN_z)
x = GN_x(i);y = GN_y(j);z = GN_z(k);
dNx = 1/8*[-(1-y)*(1-z), (1-y)*(1-z), (1+y)*(1-z), -(1+y)*(1-z), -(1-y)*(1+z), (1-y)*(1+z), (1+y)*(1+z) -(1+y)*(1+z)];
dNy = 1/8*[-(1-x)*(1-z) -(1+x)*(1-z), (1+x)*(1-z), (1-x)*(1-z), -(1-x)*(1+z), -(1+x)*(1+z), (1+x)*(1+z), (1-x)*(1+z)];
dNz = 1/8*[-(1-x)*(1-y) -(1+x)*(1-y) -(1+x)*(1+y) -(1-x)*(1+y), (1-x)*(1-y), (1+x)*(1-y), (1+x)*(1+y), (1-x)*(1+y)];
J = [dNx;dNy;dNz]*[ -a a a -a -a a a -a ;
-b -b b b -b -b b b;
-c -c -c -c c c c c]';
G = [inv(J) zeros(3) zeros(3);zeros(3) inv(J) zeros(3);zeros(3) zeros(3) inv(J)];
dN(1,1:3:24) = dNx; dN(2,1:3:24) = dNy; dN(3,1:3:24) = dNz;
dN(4,2:3:24) = dNx; dN(5,2:3:24) = dNy; dN(6,2:3:24) = dNz;
dN(7,3:3:24) = dNx; dN(8,3:3:24) = dNy; dN(9,3:3:24) = dNz;
Be = L*G*dN;
Ke = Ke + GaussWeigh(i)*GaussWeigh(j)*GaussWeigh(k)*det(J)*(Be'*DH*Be);
end
end
end
end
% === GENERATE ELEMENT STIFFNESS MATRIX ===
function [KE] = lk_H8(E0,nu)
A = [32 6 -8 6 -6 4 3 -6 -10 3 -3 -3 -4 -8;
-48 0 0 -24 24 0 0 0 12 -12 0 12 12 12];
k = 1/144*A'*[1; nu];
K1 = [k(1) k(2) k(2) k(3) k(5) k(5);
k(2) k(1) k(2) k(4) k(6) k(7);
k(2) k(2) k(1) k(4) k(7) k(6);
k(3) k(4) k(4) k(1) k(8) k(8);
k(5) k(6) k(7) k(8) k(1) k(2);
k(5) k(7) k(6) k(8) k(2) k(1)];
K2 = [k(9) k(8) k(12) k(6) k(4) k(7);
k(8) k(9) k(12) k(5) k(3) k(5);
k(10) k(10) k(13) k(7) k(4) k(6);
k(6) k(5) k(11) k(9) k(2) k(10);
k(4) k(3) k(5) k(2) k(9) k(12)
k(11) k(4) k(6) k(12) k(10) k(13)];
K3 = [k(6) k(7) k(4) k(9) k(12) k(8);
k(7) k(6) k(4) k(10) k(13) k(10);
k(5) k(5) k(3) k(8) k(12) k(9);
k(9) k(10) k(2) k(6) k(11) k(5);
k(12) k(13) k(10) k(11) k(6) k(4);
k(2) k(12) k(9) k(4) k(5) k(3)];
K4 = [k(14) k(11) k(11) k(13) k(10) k(10);
k(11) k(14) k(11) k(12) k(9) k(8);
k(11) k(11) k(14) k(12) k(8) k(9);
k(13) k(12) k(12) k(14) k(7) k(7);
k(10) k(9) k(8) k(7) k(14) k(11);
k(10) k(8) k(9) k(7) k(11) k(14)];
K5 = [k(1) k(2) k(8) k(3) k(5) k(4);
k(2) k(1) k(8) k(4) k(6) k(11);
k(8) k(8) k(1) k(5) k(11) k(6);
k(3) k(4) k(5) k(1) k(8) k(2);
k(5) k(6) k(11) k(8) k(1) k(8);
k(4) k(11) k(6) k(2) k(8) k(1)];
K6 = [k(14) k(11) k(7) k(13) k(10) k(12);
k(11) k(14) k(7) k(12) k(9) k(2);
k(7) k(7) k(14) k(10) k(2) k(9);
k(13) k(12) k(10) k(14) k(7) k(11);
k(10) k(9) k(2) k(7) k(14) k(7);
k(12) k(2) k(9) k(11) k(7) k(14)];
KE = E0/((nu+1)*(1-2*nu))*...
[ K1 K2 K3 K4;
K2' K5 K6 K3';
K3' K6 K5' K2';
K4 K3 K2 K1'];
end