-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmultimat_dif.m
184 lines (159 loc) · 6.35 KB
/
multimat_dif.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
function multimat_dif
% This function solves the multimaterial topology optimization problem
% min 1/2 \|y-yd\|^2 +
% alpha/2\sum_i\int((u_i+u_{i+1})u(x)-u_iu_{i+1})\chi_{[u_i,u_{i+1}]}
% s.t. -\div(Gu \grad y) = f, u_1 <= u(x) <= u_d
% using the approach described in the paper
% "A convex analysis approach to multi-material topology optimization"
% by Christian Clason and Karl Kunisch, see
% http://math.uni-graz.at/mobis/publications/SFB-Report-2015-011.pdf.
%
% January 14, 2016 Christian Clason <christian.clason@uni-due.de>
% http://udue.de/clason
%% setup
% problem parameters
N = 128; % number of nodes per dimension
maxit = 300; % max number of Newton steps
alpha = 1e-3; % control cost parameter (L^2)
tmin = 1e-6; % minimal step length
ub = [1.5 1.75 2 2.25 2.5]'; % vector of control states
d = length(ub); % number of control states
% uniform grid
x = linspace(-1,1,N)'; % spatial grid points (uniform in x,y)
[xx,yy] = meshgrid(x); % coordinates of nodes
h = (x(2)-x(1)); % grid size
N2 = N*N;
tplot = @(n,f,s) tplot_(n,f,s,N,xx,yy);
dplot = @(n,f,s) dplot_(n,f,s,N,xx,yy,ub);
% operators (first-order finite differences)
I = speye(N2,N2);
D1 = spdiags(ones(N,1)*[-1 1]/h,0:1,N,N);
Dx = kron(D1,speye(N,N)); DxT = Dx';
Dy = kron(speye(N,N),D1); DyT = Dy';
G1 = spdiags(ones(N,1)*[1 1 1]/6,-1:1,N,N);
G1(1,1) = 2/6; G1(end,end)= 2/6;
G = kron(speye(N,N),G1)+kron(G1,speye(N,N));
% homogeneous Dirichlet conditions: modify matrix
bnod = unique(find(abs(xx)==1 | abs(yy)==1));
dd = ones(N2,1); dd(bnod) = 0; Dnod = spdiags(dd,0,N2,N2);
bc = @(A) Dnod*A*Dnod + (I-Dnod);
% precompute for convenience
D = [Dx;Dy]; DT = D';
dot = @(x,y) x(1:N2).*y(1:N2)+x(N2+1:end).*y(N2+1:end);
mat = @(x) spdiags(x,0,N2,N2);
mvc = @(x) spdiags([x;x],0,2*N2,2*N2);
% right-hand side
f = 10+0*xx(:); f(bnod)=0;
% target obtained from reference coefficient
ue = 1.5+(xx.^2+yy.^2 < 3/4).*(xx.^2+yy.^2 > 1/4).*(xx>1/10)...
+(xx.^2+yy.^2 < 3/4).*(xx.^2+yy.^2 > 1/4).*(xx<-1/10);
z = bc(DT*mvc(G*ue(:))*D)\f;
dplot(1,G*ue(:),'reference');
tplot(2,z,'target');
%% compute control
% initialize iterates
y = randn(N2,1); % state variable
p = randn(N2,1); % adjoint state variable
as = zeros(2*d-1,N2); % active sets
% continuation
gamma = 1;
while gamma > 1e-12
fprintf('\nCompute solution for gamma = %1.3e:\n',gamma);
fprintf('It\tupdate\t\tresidual\tstep size\n');
% semismooth Newton iteration
it = 1; nold = 1e99; tau = 1; tflag = '';
ga = 1+2*gamma/alpha;
while true
% update active sets
as_old = as;
q = -G*dot(D*y,D*p);
% Q_i^gamma
as(1,:) = (q < alpha/2*(ga*ub(1)+ub(2)));
for i = 2:d-1
as(i,:) = (q > alpha/2*(ub(i-1)+ga*ub(i))) & ...
(q < alpha/2*(ga*ub(i)+ub(i+1)));
end
as(d,:) = (q > alpha/2*(ub(d-1)+ga*ub(d)));
He = as(1:d,:)'*ub;
% Q_i,i+1^gamma
for i = 1:d-1
as(d+i,:) = (q >= alpha/2*(ga*ub(i)+ub(i+1))) & ...
(q <= alpha/2*(ga*ub(i+1)+ub(i)));
He = He + (q-alpha/2*(ub(i)+ub(i+1))).*as(d+i,:)'/gamma;
end
DHe = sum(as(d+1:end,:),1)'/gamma;
% gradient
Ae = bc(DT*mvc(G*He)*D);
rhs1 = Ae*p+y-z; rhs1(bnod)=0;
rhs2 = Ae*y-f; rhs2(bnod)=0;
rhs = -[rhs1;rhs2];
nr = norm(rhs);
% line search
if nr >= nold % if no decrease: backtrack (never on first iteration)
tau = tau/2;
y = y - tau*dx(1:N2);
p = p - tau*dx(1+N2:end);
if tau < tmin % accept non-monotone step
tflag = 'n';
else % bypass rest of while loop; compute new gradient
continue;
end
end
% terminate Newton?
update = nnz((as-as_old));
fprintf('%i\t%d\t\t%1.3e\t%1.3e%c\n',...
it,update,nr,tau,tflag);
if update == 0 && nr < 1e-6 % success, solution found
break;
elseif it == maxit % failure, too many iterations
break;
end
% otherwise update information, continue
it = it+1; nold = nr; tau = 1; tflag = '';
% Newton matrix
DxP = mat(Dx*p); DyP = mat(Dy*p); DxY = mat(Dx*y); DyY = mat(Dy*y);
DHE = mat(DHe);
App = bc(-(DxT*DxP+DyT*DyP)*G*DHE*G*(DxP*Dx+DyP*Dy));
Ayy = bc(-(DxT*DxY+DyT*DyY)*G*DHE*G*(DxY*Dx+DyY*Dy));
Apy = bc(-(DxT*DxY+DyT*DyY)*G*DHE*G*(DxP*Dx+DyP*Dy));
Ayp = bc(-(DxT*DxP+DyT*DyP)*G*DHE*G*(DxY*Dx+DyY*Dy));
C = [I + App, Ae + Ayp; Ae + Apy, Ayy];
% semismooth Newton step
dx = C\rhs;
y = y + dx(1:N2);
p = p + dx(1+N2:end);
end % Newton
% check convergence
if it < maxit % converged: accept iterate
u = He; % compute control from dual iterate
dplot(3,G*u,'coefficient'); % plot coefficient
regnodes = nnz(as(d+1:end,:)); % number of nodes in Q_i,i+1^gamma
fprintf('Solution has %i node(s) in regularized active sets\n',regnodes);
gamma = gamma/2; % reduce gamma
else % not converged: reject, terminate
fprintf('Iterate rejected, returning u_gamma for gamma = %1.3e\n',gamma*2);
break;
end
end % continuation
%% plot results
dplot(3,G*u,'coefficient');
tplot(4,y,'state');
red_cost = (norm(ue(:))-norm(u(:)))/norm(ue(:)); % relative reduction in control cost
red_track = norm(y(:)-z(:))/norm(z(:)); % relative reduction in tracking cost
fprintf('Relative reduction in control cost: %1.3e, tracking: %1.3e\n',red_cost,red_track);
end % main function
function tplot_(n,f,s,N,x,y)
figure(n);
surf(x,y,reshape(f,N,N));
shading interp; lighting phong; camlight headlight; alpha(0.8);
title(s); xlabel('x_1'); ylabel('x_2');
drawnow;
end % tplot_ function
function dplot_(n,f,s,N,x,y,ub)
figure(n);
pcolor(x,y,reshape(f,N,N));
shading flat;
title(s); xlabel('x_1'); ylabel('x_2');
colorbar('Limits',[ub(1),ub(end)],'Ticks',ub); caxis([ub(1) ub(end)]);
drawnow;
end % dplot_ function