function [LocalMatrix] = LocalStiffnessMatricesPeriodic(PWL,S,QuadTriangle,l,m)


V11 = ones(QuadTriangle.number,1) *...
                       ( PWL.coefficient(1,:,l) .* PWL.coefficient(1,:,m) );
                   
V12 = ones(QuadTriangle.number,1) * ...
                       ( PWL.coefficient(2,:,l) .* PWL.coefficient(2,:,m) );
                   

V1 = S .* V11 + V12 ./ S;

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

V2 = S .* ( ones(QuadTriangle.number,1) * PWL.coefficient(1,:,l) )...
       .* PWL.evaluation(:,:,m); 

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

V3 = S .* ( PWL.evaluation(:,:,l) .* PWL.evaluation(:,:,m) );

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

LocalMatrix.Diffusion = QuadTriangle.weights * V1 ; % Diff
LocalMatrix.Advection = QuadTriangle.weights * V2 ; % Adv
LocalMatrix.Mass      = QuadTriangle.weights * V3 ; % Mass
