function RectMesh = RectMesh_Enrich2(RectMesh,BndryDescMat)
%% (Enrich2)Remedying a rectangular mesh to incorporate boundary edge info, 
% since not many mesh generators are good enough to automatically generate 
% boundary edge info.  
% 
% Here we adopt a NumBndPcs*6 Boundary Descrition Matrix (BndryDescMat,BDM): 
% Assuming domains are polygonal (boundaries are line segments), there is 
% one row for each linear pieces: 
%   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, ColoState; 2012/07--2016/09

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

for ig=1:NumEgs 
  % Utilizing the info obtained from RectMesh_Enrich1.m 
  % RectMesh.edge2elem(ig,2)>0 if it is an interior edge 
  if RectMesh.edge2elem(ig,2)>0 
    continue;
  end
  % Assuming now RectMesh.edge2elem(2,jg)==0 
  x1 = RectMesh.node(RectMesh.edge(ig,1),1);
  x2 = RectMesh.node(RectMesh.edge(ig,2),1);
  y1 = RectMesh.node(RectMesh.edge(ig,1),2);
  y2 = RectMesh.node(RectMesh.edge(ig,2),2);
  % Screening 
  for k=1:size(BndryDescMat,1)  % Number of 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 (x1-X1)*(Y2-Y1)==(X2-X1)*(y1-Y1) && (x2-X1)*(Y2-Y1)==(X2-X1)*(y2-Y1)
      RectMesh.BndryEdge(ig) = k;
      break;
    end
  end
end

return;