function [PWL] = PWL_function(mesh,QuadTriangle)


% === Compute the coefficient of the piecewise linear basic functions:
% === Exp: (c^1_1,c^1_2,c^1_3) in Phi_1 = c^1_1 x + c^1_2 y + c^1_3:

Coeff_PWLfun(1:3,:,1) = [...
    mesh.points(2,mesh.triangles(2,:)) - mesh.points(2,mesh.triangles(3,:)); ...
    mesh.points(1,mesh.triangles(3,:)) - mesh.points(1,mesh.triangles(2,:)); ...
    mesh.points(1,mesh.triangles(2,:)).* mesh.points(2,mesh.triangles(3,:))- ...
    mesh.points(1,mesh.triangles(3,:)).* mesh.points(2,mesh.triangles(2,:)) ];

Coeff_PWLfun(1:3,:,2) = [...
    mesh.points(2,mesh.triangles(3,:)) - mesh.points(2,mesh.triangles(1,:)); ...
    mesh.points(1,mesh.triangles(1,:)) - mesh.points(1,mesh.triangles(3,:)); ...
    mesh.points(1,mesh.triangles(3,:)) .* mesh.points(2,mesh.triangles(1,:))-...
    mesh.points(1,mesh.triangles(1,:)) .* mesh.points(2,mesh.triangles(3,:)) ];

Coeff_PWLfun(1:3,:,3) = [...
    mesh.points(2,mesh.triangles(1,:)) - mesh.points(2,mesh.triangles(2,:)); ...
    mesh.points(1,mesh.triangles(2,:)) - mesh.points(1,mesh.triangles(1,:)); ...
    mesh.points(1,mesh.triangles(1,:)).* mesh.points(2,mesh.triangles(2,:))- ...
    mesh.points(1,mesh.triangles(2,:)).* mesh.points(2,mesh.triangles(1,:)) ];


for j = 1:3
    for l = 1:3
        
        Coeff_PWLfun(l,:,j) = Coeff_PWLfun(l,:,j)./mesh.Area;
        
    end
end


% === Evaluate the piecewise linear basic functions in mesh.points (QuadTriangle.xnodes , QuadTriangle.ynodes)

Eval_PWLfun(:,:,1) = (...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(1,:,1)) .* QuadTriangle.xnodes + ...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(2,:,1)) .* QuadTriangle.ynodes + ...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(3,:,1)) );

Eval_PWLfun(:,:,2) = (...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(1,:,2)) .* QuadTriangle.xnodes + ...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(2,:,2)) .* QuadTriangle.ynodes + ...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(3,:,2)) );

Eval_PWLfun(:,:,3) = (...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(1,:,3)) .* QuadTriangle.xnodes + ...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(2,:,3)) .* QuadTriangle.ynodes + ...
    (ones(QuadTriangle.number,1) * Coeff_PWLfun(3,:,3)));


PWL.coefficient= Coeff_PWLfun;
PWL.evaluation = Eval_PWLfun;


end