function RectMesh = RectMesh_Enrich1(RectMesh)
%% (Enrich1) Enriching a rectangular mesh with secondary info such as 
%   edge, edge2elem, elem2edge, area, LenEg 
% Approach 1: See also BahCar_CMAM_2005
% James Liu, ColoState; 2012/07--2016/09

%% Mesh info 
NumNds = size(RectMesh.node,1); 
NumEms = size(RectMesh.elem,1); 

%% Auxiliary info: node-node to element 
NodeElem = sparse(NumNds,NumNds);
for ie=1:NumEms
  LblVrtx = RectMesh.elem(ie,1:4);
  NodeElem(LblVrtx(1),LblVrtx(2)) = ie;
  NodeElem(LblVrtx(2),LblVrtx(3)) = ie;
  NodeElem(LblVrtx(3),LblVrtx(4)) = ie;
  NodeElem(LblVrtx(4),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" 
RectMesh.edge = [I,J];

%% Secondary mesh info: edge vs its 2 neighboring elements 
RectMesh.edge2elem = zeros(NumEgs,2);
for ig=1:NumEgs
  LblVrtx = RectMesh.edge(ig,1:2);
  RectMesh.edge2elem(ig,1) = NodeElem(LblVrtx(1),LblVrtx(2));
  RectMesh.edge2elem(ig,2) = NodeElem(LblVrtx(2),LblVrtx(1));
end

% Adjusting
ig = find(RectMesh.edge2elem(:,1)>RectMesh.edge2elem(:,2));
tmp = RectMesh.edge2elem(ig,1);
RectMesh.edge2elem(ig,1) = RectMesh.edge2elem(ig,2);
RectMesh.edge2elem(ig,2) = tmp;
%
ig = find(RectMesh.edge2elem(:,1)==0);
RectMesh.edge2elem(ig,1) = RectMesh.edge2elem(ig,2);
RectMesh.edge2elem(ig,2) = 0;

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

%% Geometric info on element areas and edge lengths
%
% area for all elements
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);
RectMesh.area = abs((x2-x1).*(y2-y1));
% centers for all elements 
RectMesh.EmCntr = zeros(NumEms,2);
RectMesh.EmCntr(:,1) = 0.5*(x1+x2);
RectMesh.EmCntr(:,2) = 0.5*(y1+y2);
%
% lengths for all edges
k1 = RectMesh.edge(:,1);  k2 = RectMesh.edge(:,2);
x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2);  
x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);  
RectMesh.LenEg = sqrt((x2-x1).^2+(y2-y1).^2);

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

return;