MassSpringDamper.m

%% Parallel mass-spring-damper simulation
%
%  Author : Auralius Manurung ([e-mail])
%  Note   : 


%% Clear everything
clf;
clc;
clear all;
h = plot(0, 0);
xlabel('meter');
ylabel('meter');
mText = uicontrol('style','text');

%% Parameters
row = 4;
col = 7;
stiffness = -60;    % N/m
damping =  -0.5;    % Ns/m
mass = 0.1;         % Kg
ts = 0.005;         % Seconds
iIteration = 10000;

% Each node has 8 identical links and each link contributes to damping
damping_ = damping * 8; 

%% Build the nodes 
for c = 1: col
    for r = 1 : row
        node(r,c).initalPos = [(c - 1) / 100 (r - 1) / 100 ]; % 2cm step
        node(r,c).pos = node(r,c).initalPos;
        node(r,c).acc = [0 0];
        node(r,c).vel = [0 0];
        
        % The last row is fixed
        if (r == 1)
            node(r,c).isFixed = 1;
        else
            node(r,c).isFixed = 0;
        end
        
        % Manipulated node, position is imposed, force is generated
        node(r,c).isPositionImposed = 0; 
    end
end

% This is manpulated node 
node(row, 4).isPositionImposed = 1; 

%% Graphic thingy 
index = 1;
for c = 1 : col
    for r = 1 : row
        canvas(index,:) = node(r, c).pos;
        index = index + 1;
    end
end

canvas_min = min(canvas) - 0.03;
canvas_max = max(canvas) + 0.03;

xlim([canvas_min(1) canvas_max(1)])
ylim([canvas_min(2) canvas_max(2)])

%% Recursive update

tic % Start the time measurement

dir = 1;
for t = 1 : iIteration
    
    % Simulation time
    set(mText,'String', ts * t);
    
    % Imposed node moves up and down ever 300 iterations
    if rem(t, 300) == 0
        dir =~ dir; 
    end
    
    if (dir == 0)
        node(row, 4).pos = node(row, 4).pos + [0 0.01 * ts]; % Up
    else
        node(row, 4).pos = node(row, 4).pos - [0 0.01 * ts]; % Down
    end
    
    % Force update
    % Calculate force on each node
    for r = 1 : row
        nextRow = r + 1;
        prevRow = r - 1;
        
        for c = 1 : col
            nextCol = c + 1;
            prevCol = c - 1;
            
            f1 = [0 0];
            f2 = [0 0];
            f3 = [0 0];
            f4 = [0 0];
            f5 = [0 0];
            f6 = [0 0];
            f7 = [0 0];
            f8 = [0 0];

            % Link 1
            if (r < row && c > 1)
                l0 = node(r, c).initalPos - node(nextRow, prevCol).initalPos;
                lt = node(r, c).pos - node(nextRow, prevCol).pos;
                n = norm(lt, 2);                
                f1 = stiffness * (n - norm(l0, 2)) * lt / n;
            end

            % Link 2
            if (r < row)
                l0 = node(r, c).initalPos - node(nextRow, c).initalPos;
                lt = node(r, c).pos - node(nextRow, c).pos;
                n = norm(lt, 2);
                f2 = stiffness * (n - norm(l0, 2)) * lt / n;
            end

            % Link 3
            if (c < col)
                l0 = node(r, c).initalPos - node(r, nextCol).initalPos;
                lt = node(r, c).pos - node(r, nextCol).pos;
                n = norm(lt, 2);
                f3 = stiffness * (n - norm(l0, 2)) * lt / n;
            end

            % Link 4
            if (r > 1 && c < col)
                l0 = node(r, c).initalPos - node(prevRow, nextCol).initalPos;
                lt = node(r, c).pos - node(prevRow, nextCol).pos;
                n = norm(lt, 2);
                f4 = stiffness * (n - norm(l0, 2)) * lt / n;
            end

            % Link 5
            if (r > 1)
                l0 = node(r, c).initalPos - node(prevRow, c).initalPos;
                lt = node(r, c).pos - node(prevRow, c).pos;
                n = norm(lt, 2);
                f5 = stiffness * (n - norm(l0, 2)) * lt / n;    
            end

            % Link 6
            if (c > 1)
                l0 = node(r, c).initalPos - node(r, prevCol).initalPos;                        
                lt = node(r, c).pos - node(r, prevCol).pos; 
                n = norm(lt, 2);
                f6 = stiffness * (n - norm(l0, 2)) * lt / n;                     
            end
            
            % Link 7
            if (r < row && c < col)
                l0 = node(r, c).initalPos - node(nextRow, nextCol).initalPos;                        
                lt = node(r, c).pos - node(nextRow, nextCol).pos; 
                n = norm(lt, 2);
                f7 = stiffness * (n - norm(l0, 2)) * lt / n;                     
            end
            
            % Link 8
            if (r > 1 && c > 1)
                l0 = node(r, c).initalPos - node(prevRow, prevCol).initalPos;                        
                lt = node(r, c).pos - node(prevRow, prevCol).pos; 
                n = norm(lt, 2);
                f8 = stiffness * (n - norm(l0, 2)) * lt / n;                     
            end

            node(r,c).force =  f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + ... 
                               damping_ * node(r,c).vel;

        end
    end

    % Position, velocity, and acceelleration update
    for r = 1 : row        
        for c = 1: col
            if (node(r,c).isPositionImposed ~= 1 && node(r,c).isFixed ~= 1)
            %if (node(r,c).isFixed ~= 1)
                node(r,c).acc = node(r,c).force ./ mass;
                node(r,c).vel = node(r,c).vel + node(r,c).acc .* ts;
                node(r,c).pos = node(r,c).pos + node(r,c).vel .* ts;           
            end               
        end
    end
    
    % Graphic update       
    index = 1;    
    for c = 1 : col
        % Vertical line, going down
        for r = row : -1 : 1
            canvas(index, :) = node(r, c).pos;
            index = index + 1;
        end
        
        % Zig-zag line, going up
        for r = 1 : row
            canvas(index,:) = node(r,c).pos;
            index = index + 1;
            if (c < col)
                canvas(index ,:) = node(r, c + 1).pos;
                index = index + 1;
            end                          
        end
        
    end
    
    set(h, 'XData', canvas(:,1));
    set(h, 'YData', canvas(:,2));    
    
    drawnow;
    
end

toc % Stop the time measurement