%                 Aysajan Abidin

% Function:       [A B] = stiff_mass(p,e,t)
%                 Assembles the mass and stiffness matrices for
%                 the spatial discretization in 1D with homogeneous Dirichlet 
%                 boundary conditions.
%
% Input  
%        p:       1 x n vector of node coordinates.
%              
%        t:       3 x m matrix of element node numbers t(1,el)=nn1
%                 and t(2,el)=nn2, 1 <= el <= m, pointing to 
%                 corresponding node coordinates p(nn1) and p(nn2),
%                 and subdomain reference tags t(3,el), useful
%                 if data are given by different expresssions in
%                 different parts of the domain.
%
%        e:       2 x 2 matrix of boundary node numbers and 
%                 reference tags, with node number e(1,1)=nn pointing
%                 to the left boundary node coordinate p(nn), and e(2,2)
%                 is a reference tag (see example below).       
% Output
%        A:       The stiffness matrix, an n-x-n sparse matrix.
%        B:       The mass matrix, an n-x-n sparse matrix.   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [A,B] = stiff_mass(p,e,t)

    n = size(p, 2);
    A = zeros(n, n);
    B = zeros(n, n);
    for el = 1:size(t, 2)
        [dA, dB] = ElementContributions(p, t, el);
        nn = t(1:2, el);
        A(nn, nn) = A(nn, nn) + dA;
        B(nn, nn) = B(nn, nn) + dB;
    end
    
    A(1, 1) = A(2, 2); A(end, end) = A(2,2);
    B(1, 1) = B(2, 2); B(end, end) = B(2,2);

function [dA, dB] = ElementContributions(p, t, el)

% returns the element stiffness matrix dA and element mass 
% matrix dB for element number el defined by t(1:3,el) and
% p(t(1:2,el).

  NodeCoords = p(t(1:2, el));
  x = sum(NodeCoords)/2;             % element midpoint!
  
  dx = NodeCoords(2)-NodeCoords(1);
  
  %Guass Points, which is used when computing the masss matrix of the equations with
  %variable coefficients, i.e. take the function values evaluated at the
  %gauss quadrature points. 
  xgauss=[dx/2*(-1/sqrt(3))+x; dx/2*(1/sqrt(3))+x];
  
 % Phi = [1/2;1/2];
  phi = [(sqrt(3)+1)/(2*sqrt(3)) (sqrt(3)-1)/(2*sqrt(3));...
         (sqrt(3)-1)/(2*sqrt(3)) (sqrt(3)+1)/(2*sqrt(3)) ];
  coeff = eye(2); % Constant Coefficient, i.e, 1 in our case.
  dphi = [0, -1/6;-1/6, 0];
  Dphi = [1/dx ; -1/dx];

  dA =Dphi*Dphi'*dx;  %Element Stiffness Matrix
  
  dB = phi*coeff*phi*(dx/2); %Element Mass Matrix