%                 July 3rd, 2007 - Aysajan Abidin
%
%                 November 24th 2005 - Asta Herdis Hall

% Function:       [A B] = Assemble_stiff_mass(p,e,t,TAG)
%                 Assembles the mass and stiffness matrices for
%                 the spatial discretization.
%
% Input  
%      TAG:       Represents the spatial dimension. 
%        
%        If TAG = 1, then
%        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).       
%
%        If TAG = 2, then
%        p,e,t:   Information about the space mesh (p=points, e=edges,
%                 t=triangles). For further information try: help initmesh
% Output
%        A:       The stiffness matrix, an Np-x-Np sparse matrix.
%        B:       The mass matrix, an Np-x-Np sparse matrix.   

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

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

if TAG==1
    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);
  
   
else
    Np = size(p,2);   % number of points in the mesh
    Nt = size(t,2);   % number of triangles in the mesh

    A = sparse(Np,Np);  % The stiffness matrix
    B = sparse(Np,Np);  % The mass matrix

    for i = 1:Nt        % loop over the number of triangles
        triIdx = t(1:3,i);         % indices of triangle corner nodes
        triCoord = p(:,triIdx);    % xy-coordinates of nodes

        x = triCoord(1,:);  y = triCoord(2,:);
        dx = [y(2)-y(3);y(3)-y(1);y(1)-y(2)];
        dy = [x(3)-x(2);x(1)-x(3);x(2)-x(1)];
        triSize = 0.5 * abs(x(1)*dx(1) + x(2)*dx(2) + x(3)*dx(3));

        % Assemble matrices:
        B(triIdx,triIdx) = B(triIdx,triIdx) + triSize/12 *[2 1 1;1 2 1;1 1 2];
        A(triIdx,triIdx) = A(triIdx,triIdx) + (dx*dx' + dy*dy')/(4*triSize);
    end
end

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