How to complete sparse matrix?

Question

I am trying to solve an equation by a meshfree method. Shape functions (and their derivatives) have a compact support domain. So, over a given integration cell used for numerical integration, only a small number of shape functions will have nonzero values. At the same time, for each integration point belonging to a cell, the number and indexes of the nonzero shape functions may be different. Hence, I could not construct local stiffness matrices as in typical FEM codes, because the size of them are unknown. And it is not possible to construct some arrays and then use sparse command for assembling. I used sparse command to define stiffness matrix $K$, as follows:

K=sparse(NDOF,NDOF)

and then tried to find each element $K_{ij}$ during of computation.

for j=1:ngauss
    //subroutint to "find v" here
    [fi, fix, fiy]=shape(xg,yg,v);
    for i=1:length(fi)
        sub_K=fiy*fiy(i)+fix*fix(i));                 
        K(v(i),v(:))=K(v(i),v(:))+sparse(weight(j)*areaj*sub_K);
    end
end

where v is the array of indices of nonzero shape functions. But this kind of using sparse command is too slow. Are there any suggestions or any other way for that?

Answer 1

There are two simple solutions to speed up your matrix assembly:

Approach 1

As Kirill mentioned, you might want to use the following form sparse(i,j,v), where i,j,and v are vectors defining K(i(k),j(k)) = v(k). That will require you to rewrite a portion of your code.

Approach 2

Another approach would be to preallocate the sparse matrix storage. What happens, is that when you use K=sparse(NDOF,NDOF), no memory is reserved for your matrix $K$, so, every time a new element is added to the sparse matrix - reallocation happens (not entirely true, but pretty close).

You may want to use K = spalloc(NDOF,NDOF,NNZ) (See, Matlab help spalloc) . That will allocate the storage right from the start and avoid unnecessary reallocations inside the loop. However, you have to find or estimate NNZ (number of non-zero entries). If you are fine using a little bit extra memory, you can find $\max$ (upper bound — no reallocations will occur) or avarage (safe choice) #non-zero shape functions for every point and use it. Or you can find the exact number of NNZ elements by running your code without evaluation of the shape functions and assigning values to matrix elements.