function TriMesh = TriMesh_RefinedFrom_QuadriMesh(QuadriMesh,dir)
%% A triangular mesh is refined from a given quadrilateral mesh 
% dir = 0 specifies the inserted diagonals are forward
% dir~= 0 specifies the inserted diagonals are backward
% Graham Harper, James Liu, ColoState; 2012/07--2018/11

%% Mesh info of the quadrilateral mesh 
NDQ = QuadriMesh.NumNds;  % Number of nodes of given quadrilateral mesh 
NEQ = QuadriMesh.NumEms;  % Number of elements 
NGQ = QuadriMesh.NumEgs;  % Number of edges 
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);

%% Triangular mesh: Setup
TriMesh.flag = 1;
TriMesh.NumNds = NDQ;
TriMesh.NumEms = 2*NEQ;
TriMesh.NumEgs = NGQ + 1*NEQ;  % Add 1 edge per element

%% Triangle mesh: node coordinates 
TriMesh.node = QuadriMesh.node;  % nodes do not change

%% Triangle mesh: element vs nodes (3 vertices) 
% for the quadri with CW vertices 1,2,3,4 construct 2 triangles
% one with vertices 1,3,4 and another with vertices 2,3,1 if dir==0
% otherwise, once with vertices 1,2,4 and one with vertices 2,3,4
% 
TriMesh.elem = zeros(TriMesh.NumEms,3); 
% One quadrilateral is divided into 2 triangles and 
% 
% The 1st vertex of each triangle
TriMesh.elem(2*(1:NEQ),1) = k1; 
TriMesh.elem(2*(1:NEQ)-1,1) = k2;
% 
% The 2nd vertex of each triangle
if (dir==0)
  TriMesh.elem(2*(1:NEQ),2) = k3;
else
  TriMesh.elem(2*(1:NEQ),2) = k2;
end
TriMesh.elem(2*(1:NEQ)-1,2) = k3;
% 
% The 3rd vertex of each triangle
TriMesh.elem(2*(1:NEQ),3) = k4;
if (dir==0)
  TriMesh.elem(2*(1:NEQ)-1,3) = k1;
else
  TriMesh.elem(2*(1:NEQ)-1,3) = k4;
end

%% Triangle mesh: edge vs nodes
TriMesh.edge = zeros(TriMesh.NumEgs,2);
TriMesh.edge(1:NGQ,:) = QuadriMesh.edge;

%% Inserted edges: Either edge [1,3] or edge [2,4] for each quadrilateral 
if (dir==0)
  TriMesh.edge(NGQ+(1:NEQ),:) = QuadriMesh.elem(1:NEQ,[1,3]);
else
  TriMesh.edge(NGQ+(1:NEQ),:) = QuadriMesh.elem(1:NEQ,[2,4]);
end

%% Auxiliary sparse matrix: (node,node)->edge 
NumNds = TriMesh.NumNds;
II = TriMesh.edge(:,1);
JJ = TriMesh.edge(:,2);
NdNdEg = sparse(II,JJ,(1:TriMesh.NumEgs),NumNds,NumNds)...
       + sparse(JJ,II,(1:TriMesh.NumEgs),NumNds,NumNds);

%% Triangular mesh: element-vs-edges 
TriMesh.elem2edge = zeros(TriMesh.NumEms,3);
k1 = TriMesh.elem(:,1);  
k2 = TriMesh.elem(:,2);
k3 = TriMesh.elem(:,3);
for ie=1:TriMesh.NumEms
  TriMesh.elem2edge(ie,1) = NdNdEg(k2(ie),k3(ie));
  TriMesh.elem2edge(ie,2) = NdNdEg(k3(ie),k1(ie));
  TriMesh.elem2edge(ie,3) = NdNdEg(k1(ie),k2(ie));
end

%% Triangular mesh: edge-vs-elements 
TriMesh.edge2elem = zeros(TriMesh.NumEgs,2);
cnt = zeros(TriMesh.NumEgs,1);
for ie=1:TriMesh.NumEms
  for j=1:3
    k = TriMesh.elem2edge(ie,j);
    cnt(k) = cnt(k) + 1;
    TriMesh.edge2elem(k,cnt(k)) = ie;
  end
end

%% Triangle mesh: area 
k1 = TriMesh.elem(:,1);  
k2 = TriMesh.elem(:,2);  
k3 = TriMesh.elem(:,3); 
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2);
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
% TriMesh.area = abs(x1.*(y2-y3)+x2.*(y3-y1)+x3.*(y1-y2))/2;
TriMesh.area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));

return;