function T=expandGrid(T)

% T=expandGrid(T)
% 
% Input:
%     T  : a triangulation with the original four basic fields
%          T.coordinates, T.elements, T.dirichlet, T.neumann
% Output:
%     T  : an expanded version of the triangulation with some new fields

T.baryc=(T.coordinates(T.elements(:,1),:)...
         +T.coordinates(T.elements(:,2),:)...
         +T.coordinates(T.elements(:,3),:))/3;
if length(T.neumann)>0
    v=T.coordinates(T.neumann(:,2),:)...
        -T.coordinates(T.neumann(:,1),:);
    T.normal=[v(:,2) -v(:,1)];           % non-normalized normals
    T.midptNeu=0.5*T.coordinates(T.neumann(:,1),:)...
        +0.5*T.coordinates(T.neumann(:,2),:);
end
x12=T.coordinates(T.elements(:,2),1)-T.coordinates(T.elements(:,1),1);  % x(2)-x(1)
y12=T.coordinates(T.elements(:,2),2)-T.coordinates(T.elements(:,1),2);  % y(2)-y(1)
x13=T.coordinates(T.elements(:,3),1)-T.coordinates(T.elements(:,1),1);  % x(3)-x(1) 
y13=T.coordinates(T.elements(:,3),2)-T.coordinates(T.elements(:,1),2);  % y(3)-y(1)
T.detB=x12.*y13-x13.*y12;

% coefficients for stiffness matrix

T.c11=(x13.^2+y13.^2)./T.detB;
T.c22=(x12.^2+y12.^2)./T.detB;
T.c12=-(x12.*x13+y12.*y13)./T.detB;

% coefficients for gradients

T.g11=(y12-y13)./T.detB;
T.g12=y13./T.detB;
T.g13=-y12./T.detB;
T.g21=(x13-x12)./T.detB;
T.g22=-x13./T.detB;
T.g23=x12./T.detB;

return