%% JL20151014: TO BE REVISED FOR EFFICIENCY 
function PlygnMesh = PlygnMesh_Enrich1(PlygnMesh)
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/08 

%% Mesh info 
NumNds = PlygnMesh.NumNds;
NumEms = PlygnMesh.NumEms;
ElemType = PlygnMesh.ElemType;
node = PlygnMesh.node;
elem = PlygnMesh.elem;

%% For element center
EmCntr = zeros(NumEms,2);
for ie=1:NumEms
  EmCntr(ie,:) = sum(node(elem{ie}(1:ElemType(ie)),:))/ElemType(ie);
end

%% For area 
area = zeros(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  x1 = EmCntr(ie,1);
  y1 = EmCntr(ie,2);
  area(ie) = 0;
  for j=1:n
    j1 = mod(j,n) + 1;
    x2 = node(elem{ie}(j),1); 
    y2 = node(elem{ie}(j),2); 
    x3 = node(elem{ie}(j1),1); 
    y3 = node(elem{ie}(j1),2);
    ar = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
    area(ie) = area(ie) + ar;
  end
end

%% For radius/diameter 
radius = zeros(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  xc = EmCntr(ie,1);
  yc = EmCntr(ie,2);
  dist = zeros(n,1);
  for j=1:n
    xj = node(elem{ie}(j),1); 
    yj = node(elem{ie}(j),2);    
    dist(j) = sqrt((xj-xc)^2+(yj-yc)^2);
  end
  radius(ie) = max(dist);
end

%% Generating auxiliary mesh info
% For node + node -> element (actually cell) 
NodeElem = sparse(NumNds, NumNds);
for ie=1:NumEms
  n = ElemType(ie);
  LblVrtx = elem{ie};
  for m=1:n
    m1 = mod(m,n) + 1;
    NodeElem(LblVrtx(m), LblVrtx(m1)) = ie;
  end
end
% For node + node -> pseudo-edge 
[I,J] = find(triu(NodeElem+NodeElem'));
NumPseudoEgs = size(I,1);
NodeEdge = sparse(I,J, [1:NumPseudoEgs], NumNds,NumNds); 
NodeEdge = NodeEdge + NodeEdge';

%% Generating secondary mesh info 
% 
% For edge vs its 2 nodes 
% Note: starting node label < ending node label, since obtained from "triu" 
edge = [I,J];
NumEgs = size(edge,1);
% 
% For edge vs its 2 neighboring elements/cells
edge2elem = zeros(NumEgs,2);
LblVrtx = zeros(1,2);
for ig=1:NumEgs
  LblVrtx = edge(ig,1:2);
  edge2elem(ig,1) = NodeElem(LblVrtx(1),LblVrtx(2));
  edge2elem(ig,2) = NodeElem(LblVrtx(2),LblVrtx(1));
end
% Adjusting 
for ig=1:NumEgs
  if edge2elem(ig,1)>edge2elem(ig,2)
    tmp = edge2elem(ig,1);
    edge2elem(ig,1) = edge2elem(ig,2);
    edge2elem(ig,2) = tmp;
  end
end
for ig=1:NumEgs
  if edge2elem(ig,1)==0
    edge2elem(ig,1) = edge2elem(ig,2);
    edge2elem(ig,2) = 0;
  end
end
% 
% For element/cell vs its edges 
elem2edge = cell(NumEms,1);
for ie=1:NumEms
  n = ElemType(ie);
  LblVrtx = elem{ie};
  elem2edge{ie} = zeros(1,n);
  for m=1:n
    m1 = mod(m,n) + 1;
    elem2edge{ie}(m) = NodeEdge(LblVrtx(m),LblVrtx(m1));
  end
end 
% 
maxn = max(ElemType);

%% Finishing up 
PlygnMesh.NumEgs = NumEgs;
PlygnMesh.area = area;
PlygnMesh.EmCntr = EmCntr;
PlygnMesh.radius = radius;
PlygnMesh.maxn = maxn;
PlygnMesh.edge = edge;
PlygnMesh.elem2edge = elem2edge;
PlygnMesh.edge2elem = edge2elem;
PlygnMesh.flag = 2;

return;