function QuadriMeshNew = QuadriMesh_RefineReg(QuadriMeshOld)
%% Quadrilateral mesh: Refining regularly 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info of the old one 
NDO = QuadriMeshOld.NumNds;
NEO = QuadriMeshOld.NumEms;
NGO = QuadriMeshOld.NumEgs;
k1 = QuadriMeshOld.elem(:,1);  k2 = QuadriMeshOld.elem(:,2);  
k3 = QuadriMeshOld.elem(:,3);  k4 = QuadriMeshOld.elem(:,4);
x1 = QuadriMeshOld.node(k1,1);  y1 = QuadriMeshOld.node(k1,2);
x2 = QuadriMeshOld.node(k2,1);  y2 = QuadriMeshOld.node(k2,2);
x3 = QuadriMeshOld.node(k3,1);  y3 = QuadriMeshOld.node(k3,2);
x4 = QuadriMeshOld.node(k4,1);  y4 = QuadriMeshOld.node(k4,2);
QuadriMeshOld.EmCntr = [0.25*(x1+x2+x3+x4), 0.25*(y1+y2+y3+y4)];

%% New mesh info: Counts 
QuadriMeshNew.NumNds = NDO + NEO + NGO;
QuadriMeshNew.NumEms = 4*NEO;
QuadriMeshNew.NumEgs = 4*NEO + 2*NGO;

%% New mesh info: Node coordinates 
QuadriMeshNew.node = zeros(QuadriMeshNew.NumNds,2);
% Old nodes 
QuadriMeshNew.node(1:NDO,:) = QuadriMeshOld.node;
% New nodes on the old edges 
QuadriMeshNew.node(NDO+(1:NGO),:) = ...
  0.5*( QuadriMeshOld.node(QuadriMeshOld.edge(:,1),:) ...
      + QuadriMeshOld.node(QuadriMeshOld.edge(:,2),:) );
% New nodes in the quadrilaterals
QuadriMeshNew.node(NDO+NGO+(1:NEO),:) = QuadriMeshOld.EmCntr;

%% New mesh info: Element vs nodes 
QuadriMeshNew.elem = zeros(QuadriMeshNew.NumEms,4);
% The 1st sub-quadrilateral using the 1st vertex of the old quadrilateral 
QuadriMeshNew.elem(4*(0:NEO-1)+1,1) = QuadriMeshOld.elem(:,1);
QuadriMeshNew.elem(4*(0:NEO-1)+1,2) = NDO + QuadriMeshOld.elem2edge(:,1);
% QuadriMeshNew.elem(4*(0:NEO-1)+1,3) = NDO + NGO + QuadriMeshOld.elem2edge(:,4);
QuadriMeshNew.elem(4*(0:NEO-1)+1,3) = NDO + NGO + (1:NEO);
QuadriMeshNew.elem(4*(0:NEO-1)+1,4) = NDO + QuadriMeshOld.elem2edge(:,4);
% The 2nd sub-quadrilateral using the 2nd vertex of the old quadrilateral 
QuadriMeshNew.elem(4*(0:NEO-1)+2,1) = QuadriMeshOld.elem(:,2);
QuadriMeshNew.elem(4*(0:NEO-1)+2,2) = NDO + QuadriMeshOld.elem2edge(:,2);
% QuadriMeshNew.elem(4*(0:NEO-1)+2,3) = NDO + NGO + QuadriMeshOld.elem2edge(:,4);
QuadriMeshNew.elem(4*(0:NEO-1)+2,3) = NDO + NGO + (1:NEO);
QuadriMeshNew.elem(4*(0:NEO-1)+2,4) = NDO + QuadriMeshOld.elem2edge(:,1);
% The 3rd sub-quadrilateral using the 3rd vertex of the old quadrilateral 
QuadriMeshNew.elem(4*(0:NEO-1)+3,1) = QuadriMeshOld.elem(:,3);
QuadriMeshNew.elem(4*(0:NEO-1)+3,2) = NDO + QuadriMeshOld.elem2edge(:,3);
% QuadriMeshNew.elem(4*(0:NEO-1)+3,3) = NDO + NGO + QuadriMeshOld.elem2edge(:,4);
QuadriMeshNew.elem(4*(0:NEO-1)+3,3) = NDO + NGO + (1:NEO);
QuadriMeshNew.elem(4*(0:NEO-1)+3,4) = NDO + QuadriMeshOld.elem2edge(:,2);
% The 4th sub-quadrilateral using the 4th vertex of the old quadrilateral 
QuadriMeshNew.elem(4*(0:NEO-1)+4,1) = QuadriMeshOld.elem(:,4);
QuadriMeshNew.elem(4*(0:NEO-1)+4,2) = NDO + QuadriMeshOld.elem2edge(:,4);
% QuadriMeshNew.elem(4*(0:NEO-1)+4,3) = NDO + NGO + QuadriMeshOld.elem2edge(:,4);
QuadriMeshNew.elem(4*(0:NEO-1)+4,3) = NDO + NGO + (1:NEO);
QuadriMeshNew.elem(4*(0:NEO-1)+4,4) = NDO + QuadriMeshOld.elem2edge(:,3);

%% New mesh info: Edge vs nodes
QuadriMeshNew.edge = zeros(QuadriMeshNew.NumEgs,2);
% Odd child: Note "1"<"2"  
QuadriMeshNew.edge(2*(0:NGO-1)+1,1) = QuadriMeshOld.edge(:,1);
QuadriMeshNew.edge(2*(0:NGO-1)+1,2) = (1:NGO) + NDO;
% Even child: Note "1"<"2" 
QuadriMeshNew.edge(2*(1:NGO),1) = QuadriMeshOld.edge(:,2);
QuadriMeshNew.edge(2*(1:NGO),2) = (1:NGO) + NDO;
% QuadriMeshNew.edge(2*(1:NGO),1) = (1:NGO) + NDO;
% QuadriMeshNew.edge(2*(1:NGO),2) = QuadriMeshOld.edge(:,2);
% Four new edges inside each original quadrilateral 
% 1st new edge
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+1,1) = NDO + NGO + (1:NEO);
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+1,2) = NDO + QuadriMeshOld.elem2edge(:,1);
% 2nd new edge
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+2,1) = NDO + NGO + (1:NEO);
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+2,2) = NDO + QuadriMeshOld.elem2edge(:,2);
% 3rd new edge
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+3,1) = NDO + NGO + (1:NEO);
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+3,2) = NDO + QuadriMeshOld.elem2edge(:,3);
% 4th new edge
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+4,1) = NDO + NGO + (1:NEO);
QuadriMeshNew.edge(2*NGO+4*(0:(NEO-1))+4,2) = NDO + QuadriMeshOld.elem2edge(:,4);
% Adjusting "vertex1"<"vertex2" 
ig = find(QuadriMeshNew.edge(:,1)>QuadriMeshNew.edge(:,2));
tmp = QuadriMeshNew.edge(ig,1);
QuadriMeshNew.edge(ig,1) = QuadriMeshNew.edge(ig,2);
QuadriMeshNew.edge(ig,2) = tmp;
    
%% Auxiliary sparse matrix: (node,node)->edge 
% NdNdEg = sparse(QuadriMeshNew.NumNds, QuadriMeshNew.NumNds);
NumNds = QuadriMeshNew.NumNds;
II = QuadriMeshNew.edge(:,1);
JJ = QuadriMeshNew.edge(:,2);
NdNdEg = sparse(II,JJ,(1:QuadriMeshNew.NumEgs),NumNds,NumNds)...
       + sparse(JJ,II,(1:QuadriMeshNew.NumEgs),NumNds,NumNds);

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

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

%% Coefficients for the bilinear mapping 
k1 = QuadriMeshNew.elem(:,1);  k2 = QuadriMeshNew.elem(:,2);  
k3 = QuadriMeshNew.elem(:,3);  k4 = QuadriMeshNew.elem(:,4);
x1 = QuadriMeshNew.node(k1,1);  y1 = QuadriMeshNew.node(k1,2);
x2 = QuadriMeshNew.node(k2,1);  y2 = QuadriMeshNew.node(k2,2);
x3 = QuadriMeshNew.node(k3,1);  y3 = QuadriMeshNew.node(k3,2);
x4 = QuadriMeshNew.node(k4,1);  y4 = QuadriMeshNew.node(k4,2);
QuadriMeshNew.EmCntr = [0.25*(x1+x2+x3+x4), 0.25*(y1+y2+y3+y4)];
QuadriMeshNew.CofA = zeros(QuadriMeshNew.NumEms,4);       
QuadriMeshNew.CofA(:,1) = x1;
QuadriMeshNew.CofA(:,2) = x2-x1;
QuadriMeshNew.CofA(:,3) = x4-x1;
QuadriMeshNew.CofA(:,4) = (x1+x3)-(x2+x4);
QuadriMeshNew.CofB = zeros(QuadriMeshNew.NumEms,4);
QuadriMeshNew.CofB(:,1) = y1;
QuadriMeshNew.CofB(:,2) = y2-y1;
QuadriMeshNew.CofB(:,3) = y4-y1;
QuadriMeshNew.CofB(:,4) = (y1+y3)-(y2+y4);
QuadriMeshNew.area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
                      + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));

return;