function malha = acrescenta_celulas_na_malha(malha,Q,p)
%acrescenta aa malha (regular mas incompleta, que eh um subconjunto da malha
%uniforme com NC(1) x NC(2) celulas no dominio CD, LD) todas as celulas 
%cujas posicoes estao na matriz Q (nq x 2). 
%Os vetores ixgpontos dessas
%celulas conterao indices para o vetor p.

CD = malha.CD;
LD = malha.LD;
NC = malha.NC;
nq = size(Q,1);
assert(nq <= NC(1)*NC(2));
%Lembramos que ixc eh tal que ixc(ix+1,iy+1) eh o indice na malha
%da celula com posicao [ix,iy], ou zero se tal celula nao estah lah.
ixc = malha.ixc;
nc = nnz(ixc); %numero de celulas na malha
for i = 1:nq
    pos = Q(i,:);
    if ixc(pos(1)+1, pos(2)+1) == 0
        cel = constroi_celula(pos,NC);
        cel.ixgpontos = [];
        nc = nc + 1;
        malha.celulas(nc) = cel;
        ixc(pos(1)+1,pos(2)+1) = nc;
    end
end
assert(nc >= nq)
malha.ixc = ixc;

%distribue os pontos pelas celulas
np = size(p,1); 
assert(size(p,2) == 2)

LD_c = LD ./ NC;  %tamanho da celula
MD = CD - 0.5*LD;  %canto inferior do dominio

for i = 1:np
    xi = p(i,1);
    yi = p(i,2);
    assert((xi > MD(1)) && (xi < MD(1)+LD(1)))
    assert((yi > MD(2)) && (yi < MD(2)+LD(2)))
    posi = floor(([xi, yi] - MD) ./ LD_c);
    assert(posi(1) >= 0)
    assert(posi(1) < NC(1))
    assert(posi(2) >= 0)
    assert(posi(2) < NC(2))
    icel = ixc(posi(1)+1,posi(2)+1);
    if icel ~= 0
        %ponto estah na celula icel
        assert((icel >= 1) && (icel <= nc))
        celi = malha.celulas(icel);
        npi = size(celi.ixgpontos,2);
        celi.ixgpontos(npi+1) = i;
        malha.celulas(icel) = celi;
    end
end

