function [K,M,F]=lassema(p,t,c,a,f)
%
% function [K,M,F]=lassema(p,t,c,a,f);   
% Assembles the stiffness matrix K, mass matrix M and RHS F
% wrt the mesh defined by p and t and given functions c,a and f.

nt=size(t,2);         % Nr. of triangles
np=size(p,2);         % Nr. of points
A=pdetrg(p,t);        % areas of all triangles

%%% functions c, a and f
%%%%%%%%%%%%%%%%%%%%%%%%

% assuming the functions c, a and f are constant in each triangle
% the values calculated in the centres of mass can be used

% barycentres of triangles
tri_ind=t(1:3,:); tri_ind=tri_ind(:);
x=p(1,tri_ind)'; y=p(2,tri_ind)';         
x3=reshape(x,3,nt); y3=reshape(y,3,nt);
x=sum(x3)/3; y=sum(y3)/3;                

% functions at the barycentres
cc=eval(c);
if size(cc,2) ==1       % if c is constant: c='1'
  cc=cc*ones(1,nt);
end
aa=eval(a);
if size(aa,2) ==1
  aa=aa*ones(1,nt);
end
ff=eval(f);
if size(ff,2) ==1
  ff=ff*ones(1,nt);
end

%%% matrixes K and M

K=sparse(np,np);      % global stiffness matrix
M=sparse(np,np);      % global mass matrix
F=sparse(np,1);

for n=1:nt   % The big loop running through each triangle t.
  
  tri=t(1:3,n);  % The local current triangle, indices
  pt=p(:,tri);   %                             coordinates
  
%%% stiffness matrix K

  gx=pt(1,:); gy=pt(2,:);                    % global coords
  X1=gx(2)-gx(3);Y1=gy(2)-gy(3);             % local coords
  X2=gx(3)-gx(1);Y2=gy(3)-gy(1);
  X3=gx(1)-gx(2);Y3=gy(1)-gy(2);

  Kt=[X1*X1+Y1*Y1,X1*X2+Y1*Y2,X1*X3+Y1*Y3;   % local matrix
    X2*X1+Y2*Y1,X2*X2+Y2*Y2,X2*X3+Y2*Y3;
    X3*X1+Y3*Y1,X3*X2+Y3*Y2,X3*X3+Y3*Y3];
  
  Kt=Kt/(4*A(n))*cc(n);      
  
  % plasing the local matrix Kt into the global matrix K
  ii=tri(:,[1 1 1]);  % One way to repeat the column 3 times
  jj=ii';
       
  % Now comes the main point, the update step of K
  K=K+sparse(ii(:),jj(:),Kt(:),np,np);

%%% mass matrix M
  if any(aa)
    Mt = aa(n)*A(n)/12*(ones(3,3)+eye(3,3)); 
    M=M+sparse(ii(:),jj(:),Mt(:),np,np);
  end

%%% RHS vector F
  if any(ff)
    ft = ff(n)*A(n)/3*ones(3,1);
    F(tri)=F(tri)+ft;
  end
  
end  % end of the triangle loop



%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-05-12 09:55:21 EEST

