function [sln] = LinElas_CG_QuadriEQ1_AsmSlv(...
  EqnBC,BndryCondType,QuadriMesh,GAUSSQUAD)
%% LinElas.: CG: Enriched Lagrangian Q1 on quadri.mesh (w/ red.integration)
% Assembling & Solving 
% James Liu, ColoState; 2012/07--2020/08 

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

%% Mesh info 
NumNds = QuadriMesh.NumNds;
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);
% JL20171223: TO BE REVISED FOR EFFICIENCY !
DirichletNodeFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgs
  j1 = QuadriMesh.edge(DirichletEdge(ig),1);
  j2 = QuadriMesh.edge(DirichletEdge(ig),2);
  DirichletNodeFlag(j1) = 1;
  DirichletNodeFlag(j2) = 1;
end
DirichletNode = find(DirichletNodeFlag);

%% Setup for the discrete linear system 
DOFs = 2*NumNds + NumEgs;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);
% sln = zeros(DOFs,1);

%% Computing element matrices 
% Element matrix: strain-strain 
EltMatStrnStrn = zeros(NumEms,12,12);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  StrnPhi = zeros(NumEms,12,2,2);
  StrnPhi(:,:,1,1) = BFGrad(:,:,1,1);
  StrnPhi(:,:,2,2) = BFGrad(:,:,2,2);
  StrnPhi(:,:,1,2) = 0.5*(BFGrad(:,:,1,2)+BFGrad(:,:,2,1));
  StrnPhi(:,:,2,1) = StrnPhi(:,:,1,2);
  for i=1:12
    for j=1:12
      EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j)...
        + ( StrnPhi(:,i,1,1).*StrnPhi(:,j,1,1)...
          + StrnPhi(:,i,1,2).*StrnPhi(:,j,1,2)...
          + StrnPhi(:,i,2,1).*StrnPhi(:,j,2,1)...
          + StrnPhi(:,i,2,2).*StrnPhi(:,j,2,2)) .* jac * GAUSSQUAD.RECT(k,3);
    end
  end
end
% Element matrix: div(avg)-div(avg) 
EltMatDivDiv = zeros(NumEms,12,12);
EltDivAvg = zeros(NumEms,12);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  for j=1:12
    EltDivAvg(:,j) = EltDivAvg(:,j) + BFDiv(:,j) .*jac./area * GAUSSQUAD.RECT(k,3);
  end
end
for i=1:12
  for j=1:12
    EltMatDivDiv(:,i,j) = EltDivAvg(:,i) .* EltDivAvg(:,j) .* area;
  end
end
% Altogether 
% EltMat = zeros(NumEms,12,12);
EltMat = (2*mu)*EltMatStrnStrn + lambda*EltMatDivDiv;
% Modification with signs for edge normal vectors 
for i=1:4
  for j=1:12
    EltMat(:,8+i,j) = EltMat(:,8+i,j) .* QuadriMesh.EmEgSign(:,i);
  end
end
for i=1:12
  for j=1:4
    EltMat(:,i,8+j) = EltMat(:,i,8+j) .* QuadriMesh.EmEgSign(:,j);
  end
end

%% JL20180224: THiS IS RISKY!!!
%% Assembling GlbMat 
GlbMat = sparse(DOFs,DOFs);
pos = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  pos(:,2*j-1) = 2*k - 1;
  pos(:,2*j  ) = 2*k;
  k = QuadriMesh.elem2edge(:,j);
  pos(:,8+j) = 2*NumNds + k;
end
% Now assembling
for i=1:12 
  II = pos(:,i);
  for j=1:12 
    JJ = pos(:,j);
    GlbMat = GlbMat + sparse(II,JJ, EltMat(:,i,j), DOFs,DOFs);
  end
end

%% Assembling GlbRHS 
GlbRHS = zeros(DOFs,1);
% Computing elementwise contributions of 12 basis functions 
EltRHS = zeros(NumEms,12);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  fval = EqnBC.fxnf(qp);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  for j=1:12
    EltRHS(:,j) = EltRHS(:,j) + GAUSSQUAD.RECT(k,3)...
      * jac .*( BFVal(:,j,1).*fval(:,1) + BFVal(:,j,2).*fval(:,2) );
  end  
end
% Modification for 3 edge-based basis functions 
for j=1:4
  EltRHS(:,8+j) = EltRHS(:,8+j) .* QuadriMesh.EmEgSign(:,j);
end
% Assembling contributions of elementwise 12 basis functions 
for ie=1:NumEms
  for j=1:4  % Treating vertices and edges in one loop 
    k = QuadriMesh.elem(ie,j);
    GlbRHS(2*k-1) = GlbRHS(2*k-1) + EltRHS(ie,2*j-1);
    GlbRHS(2*k  ) = GlbRHS(2*k  ) + EltRHS(ie,2*j  );
    k = QuadriMesh.elem2edge(ie,j);
    GlbRHS(2*NumNds+k) = GlbRHS(2*NumNds+k) + EltRHS(ie,8+j);
  end
end

%% JL20171231: TO BE FINISHED/REVISED 
%% Incorporating boundary conditions: Neumann as natural 
if NumNeumannEgs>0
  % For end-nodes 
  NeumannBC = zeros(NumNeumannEgs,2,2);
  j1 = QuadriMesh.edge(NeumannEdge,1);
  j2 = QuadriMesh.edge(NeumannEdge,2);
  x1 = QuadriMesh.node(j1,1);  y1 = QuadriMesh.node(j1,2);
  x2 = QuadriMesh.node(j2,1);  y2 = QuadriMesh.node(j2,2);
  LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    tNval = EqnBC.fxntN(qp);
    NeumannBC(:,1,1) = NeumannBC(:,1,1) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,1) * tNval(:,1) .* LenEgNeumann;
    NeumannBC(:,1,2) = NeumannBC(:,1,2) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,1) * tNval(:,2) .* LenEgNeumann;
    NeumannBC(:,2,1) = NeumannBC(:,2,1) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,2) * tNval(:,1) .* LenEgNeumann;
    NeumannBC(:,2,2) = NeumannBC(:,2,2) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,2) * tNval(:,2) .* LenEgNeumann;
  end
  for ig=1:NumNeumannEgs 
    j1 = QuadriMesh.edge(NeumannEdge(ig),1);
    j2 = QuadriMesh.edge(NeumannEdge(ig),2);
    GlbRHS(2*j1-1) = GlbRHS(2*j1-1) - NeumannBC(ig,1,1);
    GlbRHS(2*j1  ) = GlbRHS(2*j1  ) - NeumannBC(ig,1,2);
    GlbRHS(2*j2-1) = GlbRHS(2*j2-1) - NeumannBC(ig,2,1);
    GlbRHS(2*j2  ) = GlbRHS(2*j2  ) - NeumannBC(ig,2,2);
  end
  % For edges 
  NeumannBC = zeros(NumNeumannEgs,1);
  j1 = QuadriMesh.edge(NeumannEdge,1);
  j2 = QuadriMesh.edge(NeumannEdge,2);
  x1 = QuadriMesh.node(j1,1);  y1 = QuadriMesh.node(j1,2);
  x2 = QuadriMesh.node(j2,1);  y2 = QuadriMesh.node(j2,2);
  LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
  nml = QuadriMesh.EgNml(NeumannEdge,:);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    tNval = EqnBC.fxntN(qp);
    psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
    NeumannBC = NeumannBC + GAUSSQUAD.LINE(k,3)...
      * sum(tNval.*nml,2) * psi .* LenEgNeumann;
  end
  GlbRHS(2*NumNds+NeumannEdge) = GlbRHS(2*NumNds+NeumannEdge) - NeumannBC;
end

%% Incorporating boundary conditions: Dirichlet as essential 
sln = zeros(DOFs,1);
% Assuming NumDirichletEgs>0 
% For nodes 
uDval = EqnBC.fxnuD(QuadriMesh.node(DirichletNode,:));
sln(2*DirichletNode-1) = uDval(:,1);
sln(2*DirichletNode  ) = uDval(:,2);
% For edges 
sln(2*NumNds+DirichletEdge) = 0;

%% Reducing...
GlbRHS = GlbRHS - GlbMat*sln;

%% Flag 
flag = ones(DOFs,1);
flag(2*DirichletNode-1) = 0;
flag(2*DirichletNode  ) = 0;
flag(2*NumNds+DirichletEdge) = 0;
FreeNdEg = find(flag);

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

return;