function [sln] = LinElas_WG_QuadriP02P02AC02P0_AsmSlv( ...
  EqnBC,BndryCondType,QuadriMesh,GAUSSQUAD)
%% LinElas.: WG(P0^2,P0^2;AC0^2,P0): Assembling & Solving on a quadri.mesh 
% James Liu, ColoState; 2017/02--2021/02 

%% Equation info 
lambda = EqnBC.lambda;
    mu = EqnBC.mu;

%% Mesh info 
NumEms = QuadriMesh.NumEms;
NumEgs = QuadriMesh.NumEgs;
  area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Sorting out boundary edges: 1=Dirichlet, 2=Neumann 
DirichletEdge = find(BndryCondType(QuadriMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(QuadriMesh.BndryEdge+1)==2);

%% Boundary info 
NumDirichletEgs = size(DirichletEdge,1);
NumNeumannEgs = size(NeumannEdge,1);

%% Setting up the discrete linear system 
DOFs = 2*NumEms + 2*NumEgs;
GlbMat = sparse(DOFs,DOFs);
GlbRHS = zeros(DOFs,1);
sln = zeros(DOFs,1);

% %% Computing elementwise matrices: strain-div formulation 
% GM = Hdiv_QuadriAC02_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);
% GMST = Hdiv_QuadriAC02_NmlzPiolaBas_GramMatST(QuadriMesh,GAUSSQUAD);
% % CDWGB = zeros(NumEms,10,8);
% CDWGB = WG_QuadriP02P02AC02_CofNmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
% EltMatStrnStrn = zeros(NumEms,10,10);
% for i=1:10
%   for j=1:10
%     EltMatStrnStrn(:,i,j) = 0;
%     for k=1:8 
%       for l=1:8 
%         EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j) ...        
%           + CDWGB(:,i,k) .* (0.5*(GM(:,k,l)+GMST(:,k,l))) .* CDWGB(:,j,l);
%       end
%     end
%   end
% end
% % CDWDB = zeros(NumEms,10);
% CDWDB = WG_QuadriP02P02P0_CofNmlzBas_DiscWkDivBasFxn(QuadriMesh);
% EltMatDivDiv = zeros(NumEms,10,10);
% for i=1:10
%   for j=1:10
%     EltMatDivDiv(:,i,j) = CDWDB(:,i) .* CDWDB(:,j) .* area;
%   end
% end
% % Elementwise matrices altogether
% % EltMat = zeros(NumEms,10,10);
% EltMat = (2*mu) * EltMatStrnStrn + lambda * EltMatDivDiv;

%% Computing elementwise matrices: grad-div formulation 
% Note GM size NumEms*8*8, CDWGB size NumEms*10*8, CDWDB size NumEms*10 
GM = Hdiv_QuadriAC02_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);
CDWGB = WG_QuadriP02P02AC02_CofNmlzPiolaBas_DiscWkGradBasFxn( ... 
  QuadriMesh,GAUSSQUAD);
EltMatGradGrad = zeros(NumEms,10,10);
for i=1:10
  for j=1:10
    EltMatGradGrad(:,i,j) = 0;
    for k=1:8
      for l=1:8
        EltMatGradGrad(:,i,j) = EltMatGradGrad(:,i,j) ...
          + CDWGB(:,i,k) .* GM(:,k,l) .* CDWGB(:,j,l);
      end
    end
  end
end
CDWDB = WG_QuadriP02P02P0_CofNmlzBas_DiscWkDivBasFxn(QuadriMesh,GAUSSQUAD);
EltMatDivDiv = zeros(NumEms,10,10);
for i=1:10
  for j=1:10
    EltMatDivDiv(:,i,j) = CDWDB(:,i) .* CDWDB(:,j) .* area;
  end
end
% EltMat size NumEms*10*10 
EltMat = mu*EltMatGradGrad + (mu+lambda)*EltMatDivDiv;

%% Assembling GlbMat  
% Positions for assembly 
pos = zeros(NumEms,10);
pos(:,1) = 2*(0:NumEms-1) + 1;
pos(:,2) = 2*(0:NumEms-1) + 2;
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  pos(:,2+2*j-1) = 2*NumEms + 2*k-1;
  pos(:,2+2*j  ) = 2*NumEms + 2*k  ;  
end
% Now assembly 
for i=1:10 
  II = pos(:,i);
  for j=1:10 
    JJ = pos(:,j);
    GlbMat = GlbMat + sparse(II,JJ, EltMat(:,i,j), DOFs,DOFs);
  end
end

%% Assembling GlbRHS 
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
II1 = 2*(0:NumEms-1) + 1;
II2 = 2*(0:NumEms-1) + 2;
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  fval = EqnBC.fxnf(qp);
  GlbRHS(II1) = GlbRHS(II1) + GAUSSQUAD.RECT(k,3) * (jac .* fval(:,1));
  GlbRHS(II2) = GlbRHS(II2) + GAUSSQUAD.RECT(k,3) * (jac .* fval(:,2));
end

%% Incorporating boundary conditions: Neumann as natural 
tNavg = zeros(NumNeumannEgs,2);
if NumNeumannEgs>0
  x1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,1),1);
  x2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,2),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,1),2);
  y2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,2),2);
  LenEg = sqrt((x1-x2).^2+(y1-y2).^2);
  qp = zeros(NumNeumannEgs,2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp(:,1) = GAUSSQUAD.LINE(k,1)*x1 + GAUSSQUAD.LINE(k,2)*x2;
    qp(:,2) = GAUSSQUAD.LINE(k,1)*y1 + GAUSSQUAD.LINE(k,2)*y2;
    tNval = EqnBC.fxntN(qp);
    tNavg(:,1) = tNavg(:,1) + GAUSSQUAD.LINE(k,3)* tNval(:,1);
    tNavg(:,2) = tNavg(:,2) + GAUSSQUAD.LINE(k,3)* tNval(:,2);
  end
  II1 = 2*NumEms + 2*NeumannEdge - 1;
  II2 = 2*NumEms + 2*NeumannEdge;
  GlbRHS(II1) = GlbRHS(II1) - tNavg(:,1).*LenEg;  % NOTE: Negative sign 
  GlbRHS(II2) = GlbRHS(II2) - tNavg(:,2).*LenEg;  % NOTE: Negative sign 
end

%% Incorporating boundary conditions: Dirichlet as essential 
% Assuming NumDirichletEgs>0 
x1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,1),1);
x2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,2),1);
y1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,1),2);
y2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,2),2);
uDavg = zeros(NumDirichletEgs,2);
qp = zeros(NumDirichletEgs,2);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.LINE(k,1)*x1 + GAUSSQUAD.LINE(k,2)*x2;
  qp(:,2) = GAUSSQUAD.LINE(k,1)*y1 + GAUSSQUAD.LINE(k,2)*y2;
  uDval = EqnBC.fxnuD(qp);
  uDavg(:,1) = uDavg(:,1) + GAUSSQUAD.LINE(k,3)* uDval(:,1);
  uDavg(:,2) = uDavg(:,2) + GAUSSQUAD.LINE(k,3)* uDval(:,2);
end

%% For reducing... 
flag = ones(DOFs,1);
flag(2*NumEms+2*DirichletEdge-1) = zeros(NumDirichletEgs,1);
flag(2*NumEms+2*DirichletEdge  ) = zeros(NumDirichletEgs,1);
EmFreeEg = find(flag);

%% Reducing...
sln(2*NumEms+2*DirichletEdge-1) = uDavg(:,1);
sln(2*NumEms+2*DirichletEdge  ) = uDavg(:,2);
GlbRHS = GlbRHS - GlbMat*sln;

%% Solving the reduced global linear system directly 
sln(EmFreeEg) = GlbMat(EmFreeEg,EmFreeEg) \ GlbRHS(EmFreeEg);

return;