function [c,h] = retire(c,cplein,h,vide);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% enleve les elements de faible energie avec 
% densites intermediaires et lissage de h
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global hm hp b p e t a f
interm = .4;
% .4 <= interm <= .5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calcul de la solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[aux,nbt] = size(t);
[K,M,F,Q,G,H,R]=assempde(b,p,e,t,c,a,f);
u = assempde(K,M,F,Q,G,H,R);
compl = G'*u;
disp('compliance = ');disp(compl);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calcul du gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d = gradj(c,h,u);
dv = min(d,norm(d,inf)*.1);
figure(1);
pdeplot(p,e,t,'xydata',dv,...
 'xystyle','flat','mesh','off');
pause(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mise a jour de la densite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[ds,I] = sort(d);
md = round(nbt*vide);
h(I(md+1:nbt)) = ones(1,nbt-md)+interm;
h = h-interm; h = max(h,hm);
h = pdeprtni(p,t,h); h = pdeintrp(p,t,h);
h = sparse(1:nbt,1:nbt,h); c = cplein*h;
h = diag(h); h = full(h); h = h';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%