function TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat)
%% (Enrich2) Enriching a triangular mesh to incorporate boundary edge info, 
% since not many mesh generators are good enough to automatically generate 
% boundary edge info.  
% 
% Here we adopt a 6*NumBndryPcs Boundary Descrition Matrix 
%   BndryDescMat (or BDM) 
% Assuming domains are polygonal (boundaries are line segments), there is 
% one column for each linear piece: 
%   Row 1,2: The starting-point x,y-coordinates 
%   Row 3,4: The ending-point x,y-coordinates 
%   Row 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--2018/11 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;

%% 
j1 = TriMesh.elem(:,1);  j2 = TriMesh.elem(:,2);  j3 = TriMesh.elem(:,3);
x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2); 
x2 = TriMesh.node(j2,1);  y2 = TriMesh.node(j2,2);
x3 = TriMesh.node(j3,1);  y3 = TriMesh.node(j3,2);
% 
% areas for all elements 
TriMesh.area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
% 
% centers for all elements 
TriMesh.EmCntr = zeros(NumEms,2);
TriMesh.EmCntr(:,1) = (1/3)*(x1+x2+x3);
TriMesh.EmCntr(:,2) = (1/3)*(y1+y2+y3);

%% 
k1 = TriMesh.edge(:,1);  k2 = TriMesh.edge(:,2);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
TriMesh.LenEg = sqrt((x1-x2).^2+(y1-y2).^2);

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

return;