function TriMesh = TriMesh_Enrich1(TriMesh)
%% (Enrich1) Enriching a triangular mesh with secondary info such as 
%   edge, edge2elem, elem2edge, area, LenEg 
% References: 
% 1) BarCar_CMAM_2005: 
%    C.Bahriawati and C.Carstensen, "Three Matlab Implementations of the
%    Lowest-Order Raviart-Thomas MFEM with a Posteriori Error Control",
%    Comput. Meth. Appl. Math, 5(2005), pp.333--361.
% 2) L.Chen, iFEM, online documents 
% James Liu, ColoState; 2012/07--2016/09

%% Approach 1): See also BahCar_CMAM_2005  
NumNds = size(TriMesh.node,1); 
NumEms = size(TriMesh.elem,1); 

%% Auxiliary info: node-node to element 
% LblVrtx = zeros(1,3);
NodeElem = sparse(NumNds,NumNds);
for ie=1:NumEms
  LblVrtx = TriMesh.elem(ie,1:3);
  NodeElem(LblVrtx(1),LblVrtx(2)) = ie;
  NodeElem(LblVrtx(2),LblVrtx(3)) = ie;
  NodeElem(LblVrtx(3),LblVrtx(1)) = ie;
end

%% Auxiliary info: node-node to edge 
[I,J] = find(triu(NodeElem+NodeElem'));
NumEgs = size(I,1);
NodeEdge = sparse(I,J,[1:NumEgs],NumNds,NumNds); 
NodeEdge = NodeEdge + NodeEdge';

%% Secondary mesh info: edge vs its 2 nodes 
% Note: start-node label < end-node label, since obtained from "triu" 
TriMesh.edge = [I,J];
TriMesh.NumEgs = size(TriMesh.edge,1);

%% Secondary mesh info: edge vs its 2 neighboring elements 
% LblVrtx = zeros(1,2);
TriMesh.edge2elem = zeros(NumEgs,2);
for ig=1:NumEgs
  LblVrtx = TriMesh.edge(ig,1:2);
  TriMesh.edge2elem(ig,1) = NodeElem(LblVrtx(1),LblVrtx(2));
  TriMesh.edge2elem(ig,2) = NodeElem(LblVrtx(2),LblVrtx(1));
end
% Adjusting 
% for ig=1:NumEgs
%   if TriMesh.edge2elem(ig,1)>TriMesh.edge2elem(ig,2)
%     tmp = TriMesh.edge2elem(ig,1);
%     TriMesh.edge2elem(ig,1) = TriMesh.edge2elem(ig,2);
%     TriMesh.edge2elem(ig,2) = tmp;
%   end
% end
ig = find(TriMesh.edge2elem(:,1)>TriMesh.edge2elem(:,2));
tmp = TriMesh.edge2elem(ig,1);
TriMesh.edge2elem(ig,1) = TriMesh.edge2elem(ig,2);
TriMesh.edge2elem(ig,2) = tmp;
% for ig=1:NumEgs
%   if TriMesh.edge2elem(ig,1)==0
%     TriMesh.edge2elem(ig,1) = TriMesh.edge2elem(ig,2);
%     TriMesh.edge2elem(ig,2) = 0;
%   end
% end
ig = find(TriMesh.edge2elem(:,1)==0);
TriMesh.edge2elem(ig,1) = TriMesh.edge2elem(ig,2);
TriMesh.edge2elem(ig,2) = 0;

%% Secondary mesh info: element vs its 3 edges 
% LblVrtx = zeros(1,3);
TriMesh.elem2edge = zeros(NumEms,3);
for ie=1:NumEms
  LblVrtx = TriMesh.elem(ie,1:3);
  TriMesh.elem2edge(ie,1) = NodeEdge(LblVrtx(2),LblVrtx(3));
  TriMesh.elem2edge(ie,2) = NodeEdge(LblVrtx(3),LblVrtx(1));
  TriMesh.elem2edge(ie,3) = NodeEdge(LblVrtx(1),LblVrtx(2));
end

% %% Approach 2): See also Chen_iFEM_2009 
% 
% NE = TriMesh.NumEms;
% 
% AllEdge = sort([TriMesh.elem(:,[2,3]); TriMesh.elem(:,[3,1]);...
%   TriMesh.elem(:,[1,2])],2);
% [TriMesh.edge, i2, j] = unique(AllEdge, 'rows');
% TriMesh.elem2edge = reshape(j,NE,3);
% 
% i1(j(3*NE:-1:1)) = 3*NE:-1:1;  i1 = i1';
% k1 = ceil(i1/NE);  t1 = i1 - NE*(k1-1);
% k2 = ceil(i2/NE);  t2 = i2 - NE*(k2-1);
% 
% ix = (i1~=i2);
% 
% % neighbor = accumarray([[t1(ix),k1(ix)];[t2,k2]],[t2(ix);t1],[NT 3]);
% 
% BndryEdge = TriMesh.edge((i1==i2),:);
% 
% TriMesh.edge2elem = [t1,t2,k1,k2];

%% Geometric info on element areas, centers, and edge lengths 
% 
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
% 
% areas for all elements 
TriMesh.area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
% 
% centers for all elements 
TriMesh.EmCntr = zeros(NumEms,2);
TriMesh.EmCntr(:,1) = (1/3)*(x1+x2+x3);
TriMesh.EmCntr(:,2) = (1/3)*(y1+y2+y3);
% 
% lengths for all edges
k1 = TriMesh.edge(:,1);  k2 = TriMesh.edge(:,2);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2);  
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);  
TriMesh.LenEg = sqrt((x2-x1).^2+(y2-y1).^2);

%% Finishing enrichment of secondary mesh info 
TriMesh.flag = 2;

return;