-
Notifications
You must be signed in to change notification settings - Fork 0
/
totalDesignEx.m
249 lines (225 loc) · 8.71 KB
/
totalDesignEx.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
function [totalMat]=totalDesignEx(xdens,xcos,xsin,xcub,nelx,nely,totvolfrac,posttreat,problem)
% totvolfrac : total volume fraction
% nelx : number of cells in horizontal direction
% nely : number of cells in vertical direction
% posttreat : 1 if post treatment (better design); 0 without (faster)
microx=100; % number of elements in horizontal direction per cell
microy=100; % number of elements in vertical direction per cell
load('ExDBStructure.mat'); % microstructure database
totalMat=zeros(nely*microy,nelx*microx);
%RETREIVE THE MICROSTRUCTURE OF EACH CELL IN THE DATABASE
for i = 1:nely
for j = 1:nelx
xdensel=xdens(i,j);
xcosel=2*xcos(i,j)-1;
xsinel=2*xsin(i,j)-1;
xcubel=xcub(i,j);
xorel=atan(xsinel/xcosel)/pi; %xorel=-1 =>angle=-pi; xorel=1 =>angle=pi
%put orientation angle in [0,pi]
if xorel<0
xorel=xorel+1;
end %now xorel=0 =>angle=0; xorel=1 =>angle=pi
%take closest microstructure in the database
dev_dens=zeros(1,34);
dev_or=zeros(1,125);
dev_cub=zeros(1,63);
for ii=1:34
dev_dens(ii)=abs(parMat(1,125*63*(ii-1)+1)-xdensel);
end
idxdens=find(dev_dens==min(dev_dens));
idxdens=(idxdens-1)*125*63+1;
for ii=1:125
dev_or(ii)=abs(parMat(2,idxdens(1,1)+63*(ii-1))-xorel);
end
idxor=find(dev_or==min(dev_or));
idxor=(idxor-1)*63;
for ii=1:63
dev_cub(ii)=abs(parMat(3,idxdens(1,1)+idxor(1,1)+ii-1)-xcubel);
end
idxcub=find(dev_cub==min(dev_cub));
idxcub=idxcub-1;
idxtot=idxdens+idxor+idxcub;
localCell=Cell_4p(microx,microy,parMat(4:end,idxtot));
%insert the cell in its place in the total design matrix
totalMat(microy*(i-1)+1:microy*i,microx*(j-1)+1:microx*j)=localCell;
end
end
if posttreat==1
totalMat=posttreatment(totalMat,nelx*microx,nely*microy,totvolfrac,problem);
end
end
function totalMat01=posttreatment(totalMat,nelxm,nelym,totvolfrac,problem)
rmin=4; %post treatment filter radius
unusedRatio=1.1; %when the stress in an element is lower than the mean stress divided by that factor, it is considered as unused. % 1.1
unusedRatio2=2; %different factor used in a secound "stress filter"
unusedRatio3=20; %different factor used in a third "stress filter"
%nelxm : total number of elements in horizontal direction
%nelym : total number of elements in vertical direction
switch problem
case 'MBB'
% USER-DEFINED ACTIVE ELEMENTS
activeelts=ones(nelxm*nelym,1);
case 'Canti'
% USER-DEFINED ACTIVE ELEMENTS
activeelts=ones(nelxm*nelym,1);
case 'Lshape'
% USER-DEFINED ACTIVE ELEMENTS
emptyelts=(nelxm/2)*(nelym)+1:(nelxm)*(nelym);
emptyelts=reshape(emptyelts, nelym,nelxm/2);
emptyelts=emptyelts(1:nelym/2,:);
emptyelts=emptyelts(:);
activeelts=ones(nelxm*nelym,1);
activeelts(emptyelts)=0;
totvolfrac=totvolfrac*3/4;
end
%%POST TREATMENT
%REMOVE "UNUSED" ELEMENTS
fprintf('Post-treatment : step 1 \n');
%group cells 4 by 4 to evaluate the stress on a coarser grid
topleftmat= totalMat(1:2:end,1:2:end);
toprightmat= totalMat(1:2:end,2:2:end);
botleftmat= totalMat(2:2:end,1:2:end);
botrightmat= totalMat(2:2:end,2:2:end);
halftotmat=(topleftmat+toprightmat+botleftmat+botrightmat)./4;
[stress]=totalFEM_S(nelxm/2,nelym/2,halftotmat,problem); % evaluate stress
stressTotMat=abs(stress).*halftotmat;
meanSTM=mean(mean(stressTotMat)); %mean stress
totalMat(repelem(stressTotMat<meanSTM/unusedRatio,2,2))=0; %eliminate "unused" elements from the total design matrix
%repeat the same process with a higher "unused" threshold
topleftmat= totalMat(1:2:end,1:2:end);
toprightmat= totalMat(1:2:end,2:2:end);
botleftmat= totalMat(2:2:end,1:2:end);
botrightmat= totalMat(2:2:end,2:2:end);
halftotmat=(topleftmat+toprightmat+botleftmat+botrightmat)./4;
[stress]=totalFEM_S(nelxm/2,nelym/2,halftotmat,problem);
stressTotMat2=abs(stress).*halftotmat;
meanSTM2=mean(mean(stressTotMat2));
totalMat(repelem(stressTotMat2<meanSTM2/unusedRatio2,2,2))=0;
%repeat the same process with a higher "unused" threshold
topleftmat= totalMat(1:2:end,1:2:end);
toprightmat= totalMat(1:2:end,2:2:end);
botleftmat= totalMat(2:2:end,1:2:end);
botrightmat= totalMat(2:2:end,2:2:end);
halftotmat=(topleftmat+toprightmat+botleftmat+botrightmat)./4;
[stress]=totalFEM_S(nelxm/2,nelym/2,halftotmat,problem);
stressTotMat3=abs(stress).*halftotmat;
meanSTM3=mean(mean(stressTotMat3));
totalMat(repelem(stressTotMat3<meanSTM3/unusedRatio3,2,2))=0;
% DENSITY FILTERING
fprintf('Post-treatment : step 2 \n');
% prepare filter
iH = ones(nelxm*nelym*(2*(ceil(rmin)-1)+1)^2,1);
jH = ones(size(iH));
sH = zeros(size(iH));
k = 0;
for i1 = 1:nelxm
for j1 = 1:nelym
e1 = (i1-1)*nelym+j1;
for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelxm)
for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nelym)
e2 = (i2-1)*nelym+j2;
k = k+1;
iH(k) = e1;
jH(k) = e2;
sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2));
end
end
end
end
H = sparse(iH,jH,sH);
Hs = sum(H,2);
totalMatF=(H*totalMat(:))./Hs; %density filtering applied here
totalMatF(:)=totalMatF(:).*activeelts; %inactive elements are set back to 0
totalMatF=reshape(totalMatF,nelym,nelxm);
%filter densities to 0 or 1 while conserving volume fraction
fprintf('Post-treatment : step 3 \n');
totalMat01=totalMatF;
lowlim=0;
highlim=1;
difVF=1;
while difVF>0.0001 && highlim-lowlim > 0.000001
midlim=(lowlim+highlim)/2;
totalMat01(totalMatF>=midlim)=1;
totalMat01(totalMatF<midlim)=0;
volfrac=mean(mean(totalMat01));
difVF=abs(volfrac-totvolfrac);
if volfrac<totvolfrac
highlim=midlim;
else
lowlim=midlim;
end
end
end
function [stress]=totalFEM_S(nelx,nely,xPhys,problem)
fsum=1.0;
%% MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
penal = 3;
%% PREPARE FINITE ELEMENT ANALYSIS
A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
D=E0/(1-nu^2)*[1 nu 0;nu 1 0; 0 0 (1-nu)/2];
B=1/2*[-1 0 1 0 1 0 -1 0;0 -1 0 -1 0 1 0 1;-1 -1 -1 1 1 1 1 -1];
DB=D*B;
Cvm=[1 -0.5 0;-0.5 1 0;0 0 3];
S=DB'*Cvm*DB;
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);
edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,nelx*nely,1);
edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1);
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
switch problem
case 'MBB'
% USER-DEFINED LOAD DOFs
loadnid = 1; % Node IDs
% loadnid = nely+1;
loaddof = 2*loadnid(:) ; % DOFs
% USER-DEFINED SUPPORT FIXED DOFs
fixednid_1 = 1:(nely+1); % Node IDs
fixednid_2 = (nelx+1)*(nely+1); % Node IDs
fixeddof = [2*fixednid_1(:)-1;2*fixednid_2(:)]; % DOFs
% USER-DEFINED ACTIVE ELEMENTS
activeelts=ones(nelx*nely,1);
case 'Canti'
% USER-DEFINED LOAD DOFs
loadnid = nelx*(nely+1)+nely/2+1; % Node IDs
loaddof = 2*loadnid(:) ; % DOFs
% USER-DEFINED SUPPORT FIXED DOFs
fixednid_1 = 1:(nely+1); % Node IDs
fixednid_2 = fixednid_1; % Node IDs
fixeddof = [2*fixednid_1(:)-1;2*fixednid_2(:)]; % DOFs
% USER-DEFINED ACTIVE ELEMENTS
activeelts=ones(nelx*nely,1);
case 'Lshape'
% USER-DEFINED LOAD DOFs
loadnid = nelx*(nely+1)+nely/2+1; % Node IDs
loaddof = 2*loadnid(:) ; % DOFs
% USER-DEFINED SUPPORT FIXED DOFs
fixednid_1 = 1:(nely+1):(nelx/2)*(nely+1)+1; % Node IDs
fixednid_2 = fixednid_1; % Node IDs
fixeddof = [2*fixednid_1(:)-1;2*fixednid_2(:)]; % DOFs
% USER-DEFINED ACTIVE ELEMENTS
emptyelts=(nelx/2)*(nely)+1:(nelx)*(nely);
emptyelts=reshape(emptyelts, nely,nelx/2);
emptyelts=emptyelts(1:nely/2,:);
emptyelts=emptyelts(:);
activeelts=ones(nelx*nely,1);
activeelts(emptyelts)=0;
end
% PREPARE FINITE ELEMENT ANALYSIS
nele = nelx*nely;
ndof = 2*(nelx+1)*(nely+1);
F = sparse(loaddof,1,-fsum,ndof,1);
U = zeros(ndof,1);
freedofs = setdiff(1:ndof,fixeddof,'stable');
xPhys(:)=xPhys(:).*activeelts;
sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),64*nele,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs) = K(freedofs,freedofs)\F(freedofs);
stress = reshape(sqrt(sum((U(edofMat)*S).*U(edofMat),2)),nely,nelx); %microscopic Von Mises Stress
end