function QuadriMesh = QuadriMesh_ConcatHoriz31( ... 
  QuadriMesh1, QuadriMesh2, QuadriMesh3) 
%% Quadri.mesh concatenation 
% Assumption: 
% James Liu, ColoState; 2014/07--2020/08 

%% Preparation and simplification 
ny = QuadriMesh1.ny;
NumNds1 = QuadriMesh1.NumNds;  NumEms1 = QuadriMesh1.NumEms;  
NumNds2 = QuadriMesh2.NumNds;  NumEms2 = QuadriMesh2.NumEms;  
NumNds3 = QuadriMesh3.NumNds;  NumEms3 = QuadriMesh3.NumEms;  

%% New/global mesh: Basic info 
QuadriMesh.NumNds = NumNds1 + NumNds2 + NumNds3 - 2*(ny+1);
QuadriMesh.NumEms = ... 
  QuadriMesh1.NumEms + QuadriMesh2.NumEms + QuadriMesh3.NumEms;
QuadriMesh.NumEgs = ... 
  QuadriMesh1.NumEgs + QuadriMesh2.NumEgs + QuadriMesh3.NumEgs - 2*ny;

%% New/global mesh: Node coordinates 
QuadriMesh.node = zeros(QuadriMesh.NumNds,2);
QuadriMesh.node(1:NumNds1,:) = QuadriMesh1.node;
QuadriMesh.node(NumNds1+[1:NumNds2-ny-1],:) = QuadriMesh2.node(ny+2:end,:); 
QuadriMesh.node(NumNds1+NumNds2-(ny+1)+[1:NumNds3-ny-1],:) = ...
  QuadriMesh3.node(ny+2:end,:); 

%% New/global mesh: Element-node info 
QuadriMesh.elem = zeros(QuadriMesh.NumEms,4);
QuadriMesh.elem(1:NumEms1,:) = QuadriMesh1.elem;
QuadriMesh.elem(NumEms1+(1:NumEms2),:)= QuadriMesh2.elem + NumNds1 - (ny+1);
QuadriMesh.elem(NumEms1+NumEms2+1:end,:) ... 
  = QuadriMesh3.elem + NumNds1 + NumNds2 - 2*(ny+1);

%% New/global mesh: edge info 

%% New/global mesh: Element centers 
% QuadriMesh.EmCntr = zeros(QuadriMesh.NumEms,2);
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.EmCntr = [0.25*(x1+x2+x3+x4), 0.25*(y1+y2+y3+y4)];

return;