function QuadriMesh = Ldom_QuadriMesh_GenInit()
%% Generating an initial quadrilateral mesh for an L-shape domain: 
% Actually Gamma-shaped [-1,1]^2\([0,1]^2): with 8 nodes, 3 elements 
% James Liu, ColoState; 2012/07--2020/11  

% The initial mesh 
NumNds = 8; 
NumEms = 3;
% % For L-domain = [-1,1]^2 \ (0,1)^2 
% node = [0 0; 0 1; -1 1; -1 0; -1 -1; 0 -1; 1,-1; 1 0];
% For (Gamma-domain) L-domain = [-1,1]^2 \ (0,1)^2 
node = [0 0; 1 0; 1 1; 0 1; -1 1; -1 0; -1,-1; 0 -1];
% Only 3 elements 
elem = [1 2 3 4; 5 6 1 4; 6 7 8 1];
QuadriMesh.NumNds = NumNds;
QuadriMesh.NumEms = NumEms;
QuadriMesh.node = node;
QuadriMesh.elem = elem;
QuadriMesh.flag = 1;

%% Auxiliary info: node-node to element 
% LblVrtx = zeros(1,4);
NodeElem = sparse(NumNds,NumNds);
for ie=1:NumEms
  LblVrtx = QuadriMesh.elem(ie,1:4);
  NodeElem(LblVrtx(1),LblVrtx(2)) = ie;
  NodeElem(LblVrtx(2),LblVrtx(3)) = ie;
  NodeElem(LblVrtx(3),LblVrtx(4)) = ie;
  NodeElem(LblVrtx(4),LblVrtx(1)) = ie;
end

%% Auxiliary info: node-node to edge 
[I,J] = find(triu(NodeElem+NodeElem'));
NumEgs = size(I,1);
NodeEdge = sparse(I,J,[1:NumEgs],NumNds,NumNds); 
NodeEdge = NodeEdge + NodeEdge';

%% Secondary mesh info: edge vs its 2 nodes 
% Note: start-node label < end-node label, since obtained from "triu" 
QuadriMesh.edge = [I,J];
QuadriMesh.NumEgs = size(QuadriMesh.edge,1);

%% Secondary mesh info: edge vs its 2 neighboring elements 
% LblVrtx = zeros(1,2);
QuadriMesh.edge2elem = zeros(NumEgs,2);
for ig=1:NumEgs
  LblVrtx = QuadriMesh.edge(ig,1:2);
  QuadriMesh.edge2elem(ig,1) = NodeElem(LblVrtx(1),LblVrtx(2));
  QuadriMesh.edge2elem(ig,2) = NodeElem(LblVrtx(2),LblVrtx(1));
end
% Adjusting so that node#1label < node#2label 
ig = find(QuadriMesh.edge2elem(:,1)>QuadriMesh.edge2elem(:,2));
tmp = QuadriMesh.edge2elem(ig,1);
QuadriMesh.edge2elem(ig,1) = QuadriMesh.edge2elem(ig,2);
QuadriMesh.edge2elem(ig,2) = tmp;
% Adjusting for boundary edges 
ig = find(QuadriMesh.edge2elem(:,1)==0);
QuadriMesh.edge2elem(ig,1) = QuadriMesh.edge2elem(ig,2);
QuadriMesh.edge2elem(ig,2) = 0;

%% Secondary mesh info: element vs its 4 edges 
% LblVrtx = zeros(1,3);
QuadriMesh.elem2edge = zeros(NumEms,4);
for ie=1:NumEms
  LblVrtx = QuadriMesh.elem(ie,1:4);
  QuadriMesh.elem2edge(ie,1) = NodeEdge(LblVrtx(1),LblVrtx(2));
  QuadriMesh.elem2edge(ie,2) = NodeEdge(LblVrtx(2),LblVrtx(3));
  QuadriMesh.elem2edge(ie,3) = NodeEdge(LblVrtx(3),LblVrtx(4));
  QuadriMesh.elem2edge(ie,4) = NodeEdge(LblVrtx(4),LblVrtx(1));
end

return;