function HybridMesh = HybridMesh_Enrich2(HybridMesh,BndryDescMat)
%% (Enrich2) Enriching a hybrid 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--2016/09

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

%% JL20160808: TO BE REVISED FOR EFFECIENCY !
for ig=1:NumEgs 
  if HybridMesh.edge2elem(ig,2)>0  % it is an interior edge 
    continue;
  end
  % Assuming now QuadriMesh.edge2elem(ig,2)==0 
  x1 = HybridMesh.node(HybridMesh.edge(ig,1),1);
  y1 = HybridMesh.node(HybridMesh.edge(ig,1),2);
  x2 = HybridMesh.node(HybridMesh.edge(ig,2),1);
  y2 = HybridMesh.node(HybridMesh.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 
      HybridMesh.BndryEdge(ig) = k;
      break;
    end
  end
end

return;