function RectMesh = RectMesh_Enrich3(RectMesh,BndryDescMat)
%% (Enrich3)Enriching a rectangular 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--2016/09

%% Retrieving basic mesh info 
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;

%% Tertiary mesh info: Element centers
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);
xc = 0.5*(x1+x2);
yc = 0.5*(y1+y2);
EmCntr = [xc, yc];

%% Teetiray mesh info: Edge normal/tangential vectors 
AllEg = RectMesh.node(RectMesh.edge(:,2),:) ...
  - RectMesh.node(RectMesh.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 RectMesh.BndryEdge(ig)>0
    NmlEg(ig,:) = BndryDescMat(RectMesh.BndryEdge(ig),5:6);
  end
end

%% Edge midpoints
EgMidPt = 0.5*(RectMesh.node(RectMesh.edge(:,1),:) ...
             + RectMesh.node(RectMesh.edge(:,2),:));

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

%% For each element, four(4) edge midpoints
xm = zeros(NumEms,4);   ym = zeros(NumEms,4);
xm(:,1) = 0.5*(x1+x2);  ym(:,1) = y1;
xm(:,3) = 0.5*(x1+x2);  ym(:,3) = y2;
xm(:,2) = x2;           ym(:,2) = 0.5*(y1+y2);
xm(:,4) = x1;           ym(:,4) = 0.5*(y1+y2);
% Signs for elementwise edge normals 
DirVec1 = [xm(:,1),ym(:,1)] - [xc,yc];
DirVec2 = [xm(:,2),ym(:,2)] - [xc,yc];
DirVec3 = [xm(:,3),ym(:,3)] - [xc,yc];
DirVec4 = [xm(:,4),ym(:,4)] - [xc,yc];
sn = zeros(NumEms,4);
sn(:,1) = sign(dot(DirVec1,NmlEg(RectMesh.elem2edge(:,1),:),2));
sn(:,2) = sign(dot(DirVec2,NmlEg(RectMesh.elem2edge(:,2),:),2));
sn(:,3) = sign(dot(DirVec3,NmlEg(RectMesh.elem2edge(:,3),:),2));
sn(:,4) = sign(dot(DirVec4,NmlEg(RectMesh.elem2edge(:,4),:),2));
% Further signs, TO BE REVISED FOR EFFICIENCY 
AuxMat = sparse(NumEms,NumEgs);
for j=1:4
  AuxMat = AuxMat + sparse((1:NumEms),RectMesh.elem2edge(:,j),...
    sn(:,j),NumEms,NumEgs);
end
sns = zeros(NumEgs,2);
for ig=1:NumEgs
  sns(ig,1) = AuxMat(RectMesh.edge2elem(ig,1),ig);
  if RectMesh.edge2elem(ig,2)>0
    sns(ig,2) = AuxMat(RectMesh.edge2elem(ig,2),ig);
  end
end

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

return;