function[Diffusion,Advection,Mass] = GlobalStiffnessMatricesPeriodic(PWL,S,mesh,QuadTriangle)  

ng = mesh.numberpoints;


Mass = sparse([],[],[],ng,ng);
Diffusion = sparse([],[],[],ng,ng);
Advection  = sparse([],[],[],ng,ng);


for l = 1:3
    for m = 1:3
        
        Ig = mesh.triangles(l,:);
        Jg = mesh.triangles(m,:);
        
        [LocalElement] = LocalStiffnessMatricesPeriodic(PWL,S,QuadTriangle,l,m); % for every triangle

        LocalElement.Diffusion = LocalElement.Diffusion .* mesh.Area; 
        LocalElement.Advection = LocalElement.Advection .* mesh.Area;
        LocalElement.Mass      = LocalElement.Mass .* mesh.Area;


        Diffusion = Diffusion + sparse(Ig,Jg,LocalElement.Diffusion,ng,ng);
        
        Advection = Advection + sparse(Ig,Jg,LocalElement.Advection,ng,ng);
        
        Mass  =  Mass  + sparse(Ig,Jg,LocalElement.Mass,ng,ng);
  
    end
end


end
