function QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat)
%% (Enrich3) Enriching a quadri. mesh with tertiary mesh info such as 
%   edge unit normal/tangential vectors, WhichEdge, etc. 
% Already enriched (flag=2), BndryDescMat needed 
% James Liu, ColoState; 2012/07--2017/10

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

%% Tertiary mesh info: Edge midpoints 
EgMidPt = zeros(NumEms,4,2);
for j=1:3
  EgMidPt(:,j,:) = 0.5*(QuadriMesh.node(QuadriMesh.elem(:,j+1),:)...
                      + QuadriMesh.node(QuadriMesh.elem(:,j),:));
end
for j=4:4
  EgMidPt(:,4,:) = 0.5*(QuadriMesh.node(QuadriMesh.elem(:,1),:)...
                      + QuadriMesh.node(QuadriMesh.elem(:,4),:));
end

%% Tertiary mesh info: Edge normal/tangential vectors 
AllEg = QuadriMesh.node(QuadriMesh.edge(:,2),:) ...
  - QuadriMesh.node(QuadriMesh.edge(:,1),:);
LenEg = sqrt(AllEg(:,1).^2+AllEg(:,2).^2);
TanEg = AllEg./[LenEg,LenEg];
NmlEg = [TanEg(:,2),-TanEg(:,1)];  % Rotating clockwise 90 degree
% JL20130801: TO BE REVISED FOR EFFICIENCY 
% Correcting normal vectors 
% for ig=1:NumEgs
%   if QuadriMesh.BndryEdge(ig)>0 
%     NmlEg(ig,:) = BndryDescMat(QuadriMesh.BndryEdge(ig),5:6);  
%   end
% end
ig = find(QuadriMesh.BndryEdge>0);
NmlEg(ig,:) = BndryDescMat(QuadriMesh.BndryEdge(ig),5:6);

%% Which edge
AuxMat = sparse(NumEms,NumEgs);
for j=1:4
  AuxMat = AuxMat + sparse((1:NumEms),QuadriMesh.elem2edge(:,j),...
    j*ones(1,NumEms),NumEms,NumEgs);
end
WhichEdge = zeros(NumEgs,2);
for ig=1:NumEgs
  WhichEdge(ig,1) = AuxMat(QuadriMesh.edge2elem(ig,1),ig);
  if QuadriMesh.edge2elem(ig,2)>0
    WhichEdge(ig,2) = AuxMat(QuadriMesh.edge2elem(ig,2),ig);
  end
end

% Signs for elementwise edge normals 
DirVec1 = squeeze(EgMidPt(:,1,:)) - QuadriMesh.EmCntr;
DirVec2 = squeeze(EgMidPt(:,2,:)) - QuadriMesh.EmCntr;
DirVec3 = squeeze(EgMidPt(:,3,:)) - QuadriMesh.EmCntr;
DirVec4 = squeeze(EgMidPt(:,4,:)) - QuadriMesh.EmCntr;
sn = zeros(NumEms,4);
% sn(:,1) = sign(dot(DirVec1,NmlEg(QuadriMesh.elem2edge(:,1),:),2));
% sn(:,2) = sign(dot(DirVec2,NmlEg(QuadriMesh.elem2edge(:,2),:),2));
% sn(:,3) = sign(dot(DirVec3,NmlEg(QuadriMesh.elem2edge(:,3),:),2));
% sn(:,4) = sign(dot(DirVec4,NmlEg(QuadriMesh.elem2edge(:,4),:),2));
% Further signs, TO BE REVISED FOR EFFICIENCY 
AuxMat = sparse(NumEms,NumEgs);
for j=1:4
  AuxMat = AuxMat + sparse((1:NumEms),QuadriMesh.elem2edge(:,j),...
    sn(:,j),NumEms,NumEgs);
end
sns = zeros(NumEgs,2);
for ig=1:NumEgs
  sns(ig,1) = AuxMat(QuadriMesh.edge2elem(ig,1),ig);
  if QuadriMesh.edge2elem(ig,2)>0
    sns(ig,2) = AuxMat(QuadriMesh.edge2elem(ig,2),ig);
  end
end

%% Finishing: flag=3 for tertiary info 
QuadriMesh.EgMidPt = EgMidPt;
QuadriMesh.LenEg = LenEg;
QuadriMesh.NmlEg = NmlEg;
QuadriMesh.TanEg = TanEg;
QuadriMesh.WhichEdge = WhichEdge;
QuadriMesh.SignEmEg = sn;
QuadriMesh.SignEgEm = sns;
QuadriMesh.flag = 3;

return;