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.
91 lines
2.3 KiB
91 lines
2.3 KiB
11 months ago
|
function [DH]=homogenization3D_me(x)
|
||
|
|
Emin = 1e-9;
|
||
|
|
x = max(x,Emin);
|
||
|
|
% x: num X num X num,
|
||
|
|
% xmin <= x <= 1 density of the material
|
||
|
|
% E0 = 1;
|
||
|
|
% nu = 0.33;
|
||
|
|
% penal=1;
|
||
|
|
|
||
|
|
% 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;];
|
||
|
|
num=size(x,1);
|
||
|
|
h=1/num;
|
||
|
|
|
||
|
|
% load KEFB_20.mat KE Fe intB
|
||
|
|
% load KEFBE_40.mat KE Fe intB
|
||
|
|
% KE = symbolicKe(num,E0,nu);%========================================可预计算
|
||
|
|
[D0,Ke,Fe] = LocalSMRT(40);
|
||
|
|
KE=Ke;
|
||
|
|
intB = GaussIntB(h);
|
||
|
|
|
||
|
|
[eleidx,mesh,VE]=periodicMesh(num);
|
||
|
|
nele=size(mesh,1);
|
||
|
|
edofMat = zeros(size(mesh,1),24);
|
||
|
|
edofMat(:,1:3:24) = mesh.*3-2;
|
||
|
|
edofMat(:,2:3:24) = mesh.*3-1;
|
||
|
|
edofMat(:,3:3:24) = mesh.*3;
|
||
|
|
|
||
|
|
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
|
||
|
|
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
|
||
|
|
% sK = reshape(KE(:)*(Emin+x(:)'.^penal*(1-Emin)),24*24*nele,1);
|
||
|
|
sK = reshape(KE(:)*(x(:).'),24*24*nele,1);
|
||
|
|
K = sparse(iK,jK,sK); K = (K+K')/2;
|
||
|
|
|
||
|
|
% Fe=symbolicFe_new(num,E0,nu);%======================================可预计算
|
||
|
|
iK = reshape(kron(edofMat,ones(6,1))',24*6*nele,1);
|
||
|
|
jK = reshape(kron(repmat([1 2 3 4 5 6],[nele,1]),ones(1,24))',24*6*nele,1);
|
||
|
|
sK = reshape(Fe(:)*(x(:).'),24*6*nele,1);
|
||
|
|
F = sparse(iK,jK,sK);
|
||
|
|
|
||
|
|
nnode=nele;% because of the periodic B.C.
|
||
|
|
|
||
|
|
% tic
|
||
|
|
u = zeros(nnode*3,6);
|
||
|
|
% u(4:end,:)=K(4:end,4:end)\F(4:end,:);
|
||
|
|
% toc
|
||
|
|
for i=1:6
|
||
|
|
u(4:end,i)=pcg(K(4:end,4:end),F(4:end,i),1e-4,300);
|
||
|
|
end
|
||
|
|
|
||
|
|
% tic
|
||
|
|
% u = zeros(nnode*3,6);
|
||
|
|
% for i = 1:6
|
||
|
|
% u(4:end,i)=pcg(K(4:end,4:end),F(4:end,i));
|
||
|
|
% end
|
||
|
|
% toc
|
||
|
|
% uRE = (u-u0)./u0;
|
||
|
|
% mean(mean(abs(uRE(4:end,:))))
|
||
|
|
|
||
|
|
% u = K\F;
|
||
|
|
% intB=compIntB(h);%==================================================可预计算
|
||
|
|
|
||
|
|
x_vec=x(:);
|
||
|
|
|
||
|
|
I=eye(6);
|
||
|
|
DH=zeros(6,6);
|
||
|
|
|
||
|
|
for iele=1:nele
|
||
|
|
element=mesh(iele,:);
|
||
|
|
edof=zeros(1,24);
|
||
|
|
edof(1:3:24)=3*element(:)-2;
|
||
|
|
edof(2:3:24)=3*element(:)-1;
|
||
|
|
edof(3:3:24)=3*element(:);
|
||
|
|
uie=u(edof,:);
|
||
|
|
De = D0*x_vec(iele);
|
||
|
|
DH=DH+De*(I*h^3-intB*uie);
|
||
|
|
% 未定义函数或变量 'I'��?
|
||
|
|
end
|
||
|
|
% RE = (DH-DH0)./DH0;
|
||
|
|
% mean(mean(abs(RE)))
|
||
|
|
% Imask = [[1,1,1,0,0,0];[1,1,1,0,0,0];[1,1,1,0,0,0];[0,0,0,1,0,0];[0,0,0,0,1,0];[0,0,0,0,0,1]];
|
||
|
|
% mean(mean(abs(RE.*Imask)))
|
||
|
|
end
|
||
|
|
|
||
|
|
|
||
|
|
|