function [Clevels,cond] = asm1dprec(fce_tau,fce_q,fce_f,a,fa,xl,el,levels)

% function [Clevels,cond] = asm1dprec(fce_tau,fce_q,xl,el,levels)

l = 0;
[Al,bl] = assemble(fce_tau,fce_q,fce_f,a,fa,xl,el,[1,max(size(xl))]);
Cl = inv(Al);

for l=1:levels
    [xl,el,El] = refine(xl,el);

    Cl = El*Cl*El';

    [Al,bl] = assemble(fce_tau,fce_q,fce_f,a,fa,xl,el,[1,max(size(xl))]);
    
    Dl = diag(diag(Al)); % Jacobi smoother
    
    Cl = Cl + inv(Dl);
end

Clevels = Cl;


% for test purposes only
e = abs(eig(Clevels*Al));
cond = max(e)/min(e);
