Vector/matrix sum of forces for mass-spring systems

Hello,

I'm working with mass-spring systems. So I have a mesh with n=1584 nodes, e=8959 elements and 2*s=21705 springs (I saved the springs in both directions). I have a nodesprings matrices NS(1584x22) with all the springs connected to each node.

for example

 NS(35,:)=[371   372   373   374   375   376   377   378   379   380   381     0     0     0     0     0     0     0     0     0     0     0]

This means that springs 371,372,373,374,... are connected to node 35.

I have the forces in the springs saved in a fs (3x21705) matrix. I have the forces because of the springs on the nodes in a fsn (3x1584) matrix.

The forces acting on the nodes because of the springs is a summation of the forces of each spring connected to this node. So for the indices in each row of NS (NS(i,:)>0), I take the sum fsn(:,i)=sum(fs(:,NS(i,indicesS)),2) with indicesS=(NS(i,:)>0), but I have to loop this for i=1:n (1584 times). The loop over all the nodes is the bottleneck of my code and I use mass-spring systems because they should be a lot faster.

Could this be done in one step?

Current code:

        MSDcube.lc=((x(1,S(1,:))-x(1,S(2,:))).^2+(x(2,S(1,:))-x(2,S(2,:))).^2+(x(3,S(1,:))-x(3,S(2,:))).^2).^(1/2); % the current length of the spring
        MSDcube.Aim=(x(1:3,S(2,:))-x(1:3,S(1,:)))./MSDcube.lc; % the direction of the springs (in 2 directions)
        fs=(k.*(MSDcube.lc-l0)).*MSDcube.Aim; %the fs matrix (3x21705) => spring forces
        for i=1:n 
            indicesS=(NS(i,:)>0); % this is done to get only the positive integers for indexaction
            MSDcube.fs(:,i)=sum(fs(:,NS(i,indicesS)),2); % MSDcube.fs is the fns (3x1584) force on each node because of the springs connected to the node.
        end
        MSDcube.fn=MSDcube.fs+MSDcube.fg;

Does anybody know if this could be done in one step through vectornotations or specific formulas? For more denser meshes this code becomes very slow because of the loop over all the nodes.

Thanks in advance,

Martijn

Answer 1

[Jaynik]

Hi Martijn,

To optimize the computation of forces on nodes in the given mass-spring system, the explicit loop over each node can be avoided by using a vectorized approach with sparse matrices. This method efficiently handles matrix operations, which can significantly speed up the computation, especially for larger and denser meshes.

The sparse matrix will represent the connectivity between nodes and springs. This matrix can then be used to perform matrix multiplication, efficiently summing the forces for each node.

Following is a sample code for the same. Feel free to modify it based on your requirement:

% Flatten the NS matrix to get non-zero entries
[rowIndices, colIndices] = find(NS > 0);
springIndices = NS(sub2ind(size(NS), rowIndices, colIndices));

% Create sparse matrix directly using these indices
node_spring_connectivity = sparse(rowIndices, springIndices, 1, n, size(fs, 2));

% Multiply the connectivity matrix by the force matrix to get the forces on nodes
MSDcube.fs = fs * node_spring_connectivity';

% Add the gravitational forces (or any other forces acting on the nodes)
MSDcube.fn = MSDcube.fs + MSDcube.fg;

The node_spring_connectivity sparse matrix efficiently maps each node to its connected springs. It is constructed using the indices of non-zero elements in NS, avoiding the need for an explicit loop.

Refer to the following documentation to read more about these functions:

• sparse: https://www.mathworks.com/help/matlab/ref/sparse.html
• sub2ind: https://www.mathworks.com/help/matlab/ref/sub2ind.html
• find: https://www.mathworks.com/help/matlab/ref/find.html

Hope this helps!