%% JL20171209: TO BE FINISHED/REVISED 
function TriMesh = TriMesh_Enrich4(TriMesh,BndryDescMat)
%% (Enrich4) Enriching a triangular mesh for more edge info 
% James Liu, ColoState; 2012/07--2018/11 

%% Mesh info 
NumNds = TriMesh.NumNds;
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;

%% Computing edge length and normal vector 
EgVec = TriMesh.node(TriMesh.edge(:,2),:)...
      - TriMesh.node(TriMesh.edge(:,1),:);
EgLen = sqrt(EgVec(:,1).^2+EgVec(:,2).^2);
EgTan = EgVec./[EgLen,EgLen];
EgNml = [EgTan(:,2),-EgTan(:,1)];  % Rotating clockwise 90 degree
% Adjusting edge normal vector for boundary edges 
kg = find(TriMesh.BndryEdge>0);
EgNml(kg,:) = BndryDescMat(TriMesh.BndryEdge(kg),5:6);
% Incorporating 
TriMesh.EgLen = EgLen;
TriMesh.EgNml = EgNml;

%% Initialization and Correction
TriMesh.EgEmSign = -ones(NumEgs,2);
EgMidPt = 0.5*(TriMesh.node(TriMesh.edge(:,2),:)...
             + TriMesh.node(TriMesh.edge(:,1),:));
DirVec1 = EgMidPt - TriMesh.EmCntr(TriMesh.edge2elem(:,1),:);
dgd = dot(DirVec1,EgNml,2);
[(1:NumEgs)', dgd];
kg = find(TriMesh.BndryEdge==0 & dgd>0);  
TriMesh.EgEmSign(kg,1) =  1;  
TriMesh.EgEmSign(kg,2) = -1;
kg = find(TriMesh.BndryEdge==0 & dgd<0);
TriMesh.EgEmSign(kg,1) = -1;  
TriMesh.EgEmSign(kg,2) =  1;
kg1 = find(TriMesh.BndryEdge>0);  
TriMesh.EgEmSign(kg1,1) = 1;  
TriMesh.EgEmSign(kg1,2) = 0;

%% 
AuxMatEgEm = sparse(NumEgs,NumEms);
II = (1:NumEgs);  
JJ1 = TriMesh.edge2elem(:,1);
AuxMatEgEm = AuxMatEgEm + sparse(II,JJ1, TriMesh.EgEmSign(:,1), NumEgs,NumEms);
II = find(TriMesh.BndryEdge==0);  
JJ2 = TriMesh.edge2elem(II,2);
AuxMatEgEm = AuxMatEgEm + sparse(II,JJ2, TriMesh.EgEmSign(II,2), NumEgs,NumEms);

%% 
TriMesh.EmEgSign = zeros(NumEms,3);  % Default values 
for ie=1:NumEms
  TriMesh.EmEgSign(ie,1) = AuxMatEgEm(TriMesh.elem2edge(ie,1),ie);
  TriMesh.EmEgSign(ie,2) = AuxMatEgEm(TriMesh.elem2edge(ie,2),ie);
  TriMesh.EmEgSign(ie,3) = AuxMatEgEm(TriMesh.elem2edge(ie,3),ie);
end

return;
