%% JL20171223: TO BE FINISHED/REVISED 
function QuadriMesh = QuadriMesh_Enrich4(QuadriMesh,BndryDescMat)
%% (Enrich4) Enriching a quadrilateral mesh for more edge info 
% Using after Enrich2 
% James Liu, ColoState; 2012/07--2017/12 

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

%% Computing edge length and normal vector 
EgVec = QuadriMesh.node(QuadriMesh.edge(:,2),:)...
      - QuadriMesh.node(QuadriMesh.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(QuadriMesh.BndryEdge>0);
EgNml(kg,:) = BndryDescMat(QuadriMesh.BndryEdge(kg),5:6);
% Incorporating 
QuadriMesh.EgLen = EgLen;
QuadriMesh.EgNml = EgNml;

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

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

%% EmEgSign: element-edge signs 
QuadriMesh.EmEgSign = zeros(NumEms,4);  % Default values 
for ie=1:NumEms
  QuadriMesh.EmEgSign(ie,1) = AuxMatEgEm(QuadriMesh.elem2edge(ie,1),ie);
  QuadriMesh.EmEgSign(ie,2) = AuxMatEgEm(QuadriMesh.elem2edge(ie,2),ie);
  QuadriMesh.EmEgSign(ie,3) = AuxMatEgEm(QuadriMesh.elem2edge(ie,3),ie);
  QuadriMesh.EmEgSign(ie,4) = AuxMatEgEm(QuadriMesh.elem2edge(ie,4),ie);
end

return;