%% This Routine is for plane stress for triangles and a banded matrix

function [Lhs] = plnstr(x,y,vu,in,dof,ne,BW)
%
% this subroutine builds the linearized equations of elasticity
%    for plane stress for triangles
% It assumes Lhs(2*number nodes, bandwidth)  
%  dof is the total number of degrees of freedom = 2*number nodes
%
%  initialize the lhs matrix:  Lhs( , )
%  element are dimensioned in(3,ne)  3=triangles  ne=number of elements
%
Lhs = zeros(dof,BW);
%
%  poissons ratio vu  (nu)
%
vu2 = (1.-vu)/2.;
hb = (BW-1)/2;    diag = hb + 1;
for L = 1:ne
    for i = 1:3
        j = in(mod(i,3)+1,L);
        k = in(mod(i+1,3)+1,L);
        dx(i) = x(j)-x(k);
        dy(i) = y(j)-y(k);
    end
    area = (-dx(3)*dy(2) + dy(3)*dx(2) )/2;
    a4 = 1./(4*area);    % "/" is the slowest basic operation
    for i=1:3
        iy = 2*in(i,L);    ix = iy-1;
        for j=1:3
            jy = 2*in(j,L);    jx = jy-1;
            k = diag + jx-ix;
            Lhs(ix,k) = Lhs(ix,k) + (dy(j)*dy(i) + vu2*dx(j)*dx(i))*a4;
            k = ndiag2 + jy-ix;
            Lhs(ix,k) = Lhs(ix,k) - (vu*dx(j)*dy(i) + vu2*dy(j)*dx(i))*a4;
            k = ndiag2 + jx-iy;
            Lhs(iy,k) = Lhs(iy,k) - (vu*dy(j)*dx(i) + vu2*dx(j)*dy(i))*a4;
            k = ndiag2 + jy-iy;
            Lhs(iy,k) = Lhs(iy,k) + (dx(j)*dx(i) + vu2*dy(j)*dy(i))*a4;
        end   % j loop
    end % i loop
end % L loop

return
end

%%  This code puts the assembled coefficients into a sparse matrix storage

function [a] = plnstr(x,y,vu,in,dof,ne,a,ia,ja)
%
% this subroutine builds the linearized equations of elasticity
%    for plane stress for triangles
% It assumes ia(dof) indicates the start location in the ja() and a() arrays
%             ja(dof*number of neighbors) identifies which columns 
%               in a virtual A[N,N] matrix is being addressed
%               and sums the coefficients into a sparse vector a( )
%  dof is the total number of degrees of freedom = 2*number nodes
%
%  element are dimensioned in(3,ne)  3=triangles  ne=number of elements
%
Lhs = zeros(dof,BW);
%
%  poissons ratio vu  (nu)
%
vu2 = (1.-vu)/2.;
hb = (BW-1)/2;    diag = hb + 1;
for L = 1:ne
    for i = 1:3
        j = in(mod(i,3)+1,L);
        k = in(mod(i+1,3)+1,L);
        dx(i) = x(j)-x(k);
        dy(i) = y(j)-y(k);
    end
    area = (-dx(3)*dy(2) + dy(3)*dx(2) )/2;
    a4 = 1./(4*area);    % "/" is the slowest basic operation
    for i=1:3
        iy = ia(2*in(i,L));    
        ix = ia(2*in(i,L)-1);
        for j=1:3
            jy = 2*in(j,L);    
            jx = jy-1;
            % sum F(x) is a fxn of (u and v); this is for u displacement
            % find jx in the ja( array ) within the row of ix
            loc = find_column(ix,jx,ia,ja);
            a(loc) = a(loc) + (dy(j)*dy(i) + vu2*dx(j)*dx(i))*a4;
            % sum F(x) is a fxn of (u and v); this is for v displacement
            loc = find_column(ix,jy,ia,ja);
            a(loc) = a(loc) - (vu*dx(j)*dy(i) + vu2*dy(j)*dx(i))*a4;
            % sum F(y) is a fxn of (u and v); this is for u displacement
            loc = find_column(iy,jx,ia,ja);
            a(loc) = a(loc) - (vu*dy(j)*dx(i) + vu2*dx(j)*dy(i))*a4;
            % sum F(y) is a fxn of (u and v); this is for v displacement
            loc = find_column(iy,jy,ia,ja);
            a(loc) = a(loc) + (dx(j)*dx(i) + vu2*dy(j)*dy(i))*a4;
        end   % j loop
    end % i loop
end % L loop

return
end

function [loc] = find_column(row, col, ia, ja)
    Row_start = ia(row);
    Row_end = ia(row+1)-1;
    for i = Row_start:Row_end
        if ja(i) == col
            loc = i;
            break
        end
    end
    % would need some flag if col was not found
end