function QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat)
%% (Enrich2) Enriching a quadri. mesh to incorporate boundary edge info, 
% since few mesh generators automatically generate such info  
% 
% Here we adopt a Boundary Descrition Matrix (BDM) of size NumBndryPcs*6 
% Assuming the domain is polygonal (boundaries are line segments), 
% there is one row in BDM for each linear piece: 
%   Column 1,2: The starting-point x,y-coordinates 
%   Column 3,4: The ending-point x,y-coordinates 
%   Column 5,6: The outward unit normal vector x,y-componets 
% The BDM is problem-dependent but the code is not.
% 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2019/11 

%%  
NumEgs = QuadriMesh.NumEgs;
QuadriMesh.BndryEdge = zeros(NumEgs,1);

%% JL20160808: TO BE REVISED FOR EFFECIENCY !
for ig=1:NumEgs 
  if QuadriMesh.edge2elem(ig,2)>0  % it is an interior edge 
    continue;
  end
  % Assuming now QuadriMesh.edge2elem(ig,2)==0 
  x1 = QuadriMesh.node(QuadriMesh.edge(ig,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(ig,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(ig,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(ig,2),2);
  % Screening 
  for k=1:size(BndryDescMat,1)  % Number of pieces (segments) on the boundary 
    X1 = BndryDescMat(k,1);  Y1 = BndryDescMat(k,2);  % Starting point 
    X2 = BndryDescMat(k,3);  Y2 = BndryDescMat(k,4);  % Ending point 
    if abs((x1-X1)*(Y2-Y1)-(X2-X1)*(y1-Y1))<1e-9 && abs((x2-X1)*(Y2-Y1)-(X2-X1)*(y2-Y1))<1e-9 
      QuadriMesh.BndryEdge(ig) = k;
      break;
    end
  end
end

%% Elementwise area, center, diameter
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
QuadriMesh.area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
                       + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));
QuadriMesh.EmCntr = [0.25*(x1+x2+x3+x4), 0.25*(y1+y2+y3+y4)];
diag13 = sqrt((x1-x3).^2+(y1-y3).^2);
diag24 = sqrt((x2-x4).^2+(y2-y4).^2);
QuadriMesh.diam = max(diag13,diag24);

return;