%% JL20160612: SOME PARTS NEED DOUBLE-CHECK!!
function QuadriMesh = QuadriMesh_RefinedFrom_TriMesh(TriMesh)
%% A quadrlateral mesh refined from a given triangular mesh 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2016/06

%% Getting info of trig.mesh 
NDT = TriMesh.NumNds;  % Number of nodes of the given triangular mesh 
NET = TriMesh.NumEms;  % Number of elements
NGT = TriMesh.NumEgs;  % Number of edges
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);
xc = (x1+x2+x3)/3;
yc = (y1+y2+y3)/3;

%% Quadrilateral mesh: counts 
QuadriMesh.flag = 1;
QuadriMesh.NumNds = NDT + NET + NGT;
QuadriMesh.NumEms = 3*NET;
QuadriMesh.NumEgs = 3*NET + 2*NGT;
% QuadriMesh.node = zeros(QuadriMesh.NumNds,2);
% QuadriMesh.elem = zeros(QuadriMesh.NumEms,4);
% QuadriMesh.edge = zeros(QuadriMesh.NumEgs,2);

%% Quadrilateral mesh: node coordinates  
QuadriMesh.node = zeros(QuadriMesh.NumNds,2);
% The 1st group of nodes 
QuadriMesh.node(1:NDT,:) = TriMesh.node;  % Original nodes of tri.mesh 
% The 2nd group of nodes: center of the original triangles
QuadriMesh.node(NDT+(1:NET),1) = xc;
QuadriMesh.node(NDT+(1:NET),2) = yc;
% The 3rd group of nodes: midpoints of the original edges
QuadriMesh.node(NDT+NET+(1:NGT),:) = ...
  0.5*(TriMesh.node(TriMesh.edge(:,1),:)+TriMesh.node(TriMesh.edge(:,2),:));

%% Quadrilateral mesh: element vs nodes (its 4 vertices) 
QuadriMesh.elem = zeros(QuadriMesh.NumEms,4); 
% One triangle is divided into 3 quadrilaterals and 
% The 1st vertex of each quadrilateral = a vertex of original triangle 
QuadriMesh.elem(3*(0:NET-1)+1,1) = k1; 
QuadriMesh.elem(3*(0:NET-1)+2,1) = k2;
QuadriMesh.elem(3*(0:NET-1)+3,1) = k3;
% The 2nd vertex of each quadrilateral 
QuadriMesh.elem(3*(0:NET-1)+1,2) = (NDT+NET) + TriMesh.elem2edge(:,3);
QuadriMesh.elem(3*(0:NET-1)+2,2) = (NDT+NET) + TriMesh.elem2edge(:,1);
QuadriMesh.elem(3*(0:NET-1)+3,2) = (NDT+NET) + TriMesh.elem2edge(:,2);
% The 3rd vertex of each quadrilateral = center of the original triangle 
QuadriMesh.elem(3*(0:NET-1)+1,3) = NDT + (1:NET);
QuadriMesh.elem(3*(0:NET-1)+2,3) = NDT + (1:NET);
QuadriMesh.elem(3*(0:NET-1)+3,3) = NDT + (1:NET);
% The 4th vertex of each quadrilateral 
QuadriMesh.elem(3*(0:NET-1)+1,4) = (NDT+NET) + TriMesh.elem2edge(:,2);
QuadriMesh.elem(3*(0:NET-1)+2,4) = (NDT+NET) + TriMesh.elem2edge(:,3);
QuadriMesh.elem(3*(0:NET-1)+3,4) = (NDT+NET) + TriMesh.elem2edge(:,1);

%% Quadrilateral mesh: edge vs nodes (its 2 vertices) 
% The original edge is split into 2 child edges 
QuadriMesh.edge = zeros(QuadriMesh.NumEgs,2);
% Left child 1st vertex  = original edge 1st vertex 
QuadriMesh.edge(2*(0:NGT-1)+1,1) = TriMesh.edge(:,1);
% Right child 2nd vertex = original edge 2nd vertex 
QuadriMesh.edge(2*(1:NGT),2) = TriMesh.edge(:,2);
% Left child 2nd vertex = Right child 1st vertex = midpt of original edge 
QuadriMesh.edge(2*(0:NGT-1)+1,2) = (NDT+NET) + (1:NGT);
QuadriMesh.edge(2*(1:NGT),1)     = (NDT+NET) + (1:NGT);
% Three (3) new edges inside original triangle 
% 1st new edge opposite to the 1st original triangle vertex 
QuadriMesh.edge(2*NGT+3*(0:NET-1)+1,1) = (NDT+NET) + TriMesh.elem2edge(:,1);
QuadriMesh.edge(2*NGT+3*(0:NET-1)+1,2) = NDT + (1:NET);
% 2nd new edge opposite to the 2nd original triangle vertex 
QuadriMesh.edge(2*NGT+3*(0:NET-1)+2,1) = (NDT+NET) + TriMesh.elem2edge(:,2);
QuadriMesh.edge(2*NGT+3*(0:NET-1)+2,2) = NDT + (1:NET);
% 3rd new edge opposite to the 3rd original triangle vertex 
QuadriMesh.edge(2*NGT+3*(0:NET-1)+3,1) = (NDT+NET) + TriMesh.elem2edge(:,3);
QuadriMesh.edge(2*NGT+3*(0:NET-1)+3,2) = NDT + (1:NET);

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

%% New mesh info: Elements-vs-edges 
QuadriMesh.elem2edge = zeros(QuadriMesh.NumEms,4);
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
for ie=1:QuadriMesh.NumEms
  QuadriMesh.elem2edge(ie,1) = NdNdEg(k1(ie),k2(ie));
  QuadriMesh.elem2edge(ie,2) = NdNdEg(k2(ie),k3(ie));
  QuadriMesh.elem2edge(ie,3) = NdNdEg(k3(ie),k4(ie));
  QuadriMesh.elem2edge(ie,4) = NdNdEg(k4(ie),k1(ie));
end

%% New mesh info: Edges-vs-elements 
QuadriMesh.edge2elem = zeros(QuadriMesh.NumEgs,2);
cnt = zeros(QuadriMesh.NumEgs,1);
for ie=1:QuadriMesh.NumEms
  for j=1:4
    k = QuadriMesh.elem2edge(ie,j);
    cnt(k) = cnt(k) + 1;
    QuadriMesh.edge2elem(k,cnt(k)) = ie;
  end
end

%% Computing area 
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*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1))...
                + 0.5*abs((x3-x1).*(y4-y1)-(x4-x1).*(y3-y1));

return;