-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmultimat_pot.m
161 lines (139 loc) · 5.62 KB
/
multimat_pot.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
function multimat_pot
% 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. -\Delta y + u = f, \partial_\nu y = 0, 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-6; % control cost parameter (L^2)
tmin = 1e-6; % minimal step length
ub = [1 1.5 2 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
h2 = (x(2)-x(1))^2; % grid size (squared)
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);
% differential operator (second order finite differences for Poisson)
D2 = spdiags(ones(N,1)*[-1 2 -1]/h2,-1:1,N,N);
D2(1,1) = 1/h2; D2(end,end) = 1/h2; % modify for Neumann b.c.
A = kron(speye(N,N),D2)+kron(D2,speye(N,N));
% right-hand side
f = sin(pi*xx(:)).*sin(pi*yy(:));
% 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 = (A+spdiags(ue(:),0,N2,N2))\f;
dplot(1,ue,'reference');
tplot(2,z,'target');
%% compute control
% initialize iterates
y = zeros(N2,1); % state variable
p = zeros(N2,1); % adjoint state variable
as = zeros(2*d-1,N2); % active sets
% continuation strategy
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
% compute active sets
as_old = as;
q = -p(:).*y(:);
% 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
rhs = -[A*p+He.*p+(y-z(:)); A*y+He.*y-f];
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
Ae = A+spdiags(He-DHe.*p.*y,0,N2,N2);
C = [spdiags(1-DHe.*p.^2,0,N2,N2) Ae; ...
Ae' spdiags(-DHe.*y.^2,0,N2,N2)];
% 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, continue
u = He; % compute control from dual iterate
dplot(3,u,'control'); % plot control
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,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