function [sln,GlbMat,GlbRHS] = Darcy_WG_QuadriP1P0P02_AsmSlv(...
  EqnBC,BndryCondType,QuadriMesh,PermK,GAUSSQUAD,rho)
%% Darcy: WG(P1,P0;P0^2)Quadri (Stabilizer needed): Assembling and Solving 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% Cf. Mu,Wang,Ye, JCAM, 2015 
% JL20170212: TO BE FINISHED/REVISED 
% James Liu, ColoState; 2012/07--2018/06 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
NumEgs = QuadriMesh.NumEgs;
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);
area = 0.5*( abs((x2-x1).*(y3-y1) - (x3-x1).*(y2-y1))...
           + abs((x3-x1).*(y4-y1) - (x4-x1).*(y3-y1)));
xc = 0.25*(x1+x2+x3+x4);
yc = 0.25*(y1+y2+y3+y4);
diag13 = sqrt((x1-x3).^2+(y1-y3).^2);
diag24 = sqrt((x2-x4).^2+(y2-y4).^2);
diam = max(diag13,diag24);
dm1 = 1./diam;
CofA = zeros(NumEms,4);       CofB = zeros(NumEms,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);

%% 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);

%% Setup 
DOFs = 3*NumEms + NumEgs;
GlbMat = sparse(DOFs,DOFs);
% GlbStfMat = sparse(DOFs,DOFs);
% GlbStabMat = sparse(DOFs,DOFs);
% EltStfMat = zeros(NumEms,7,7);
% EltStabMat = zeros(NumEms,7,7);
% GlbRHS = zeros(DOFs,1);

%% Computing discrete weak gradients 
CDWGB = WG_QuadriP1P0P02_CofNatBas_DiscWkGradBasFxn(QuadriMesh);

%% Computing EltStfMat: Elementwise Stiffness Matrix  
EltStfMat = zeros(NumEms,7,7);
for i=4:7
  for j=4:7
    EltStfMat(:,i,j) = (PermK(:,1,1).*CDWGB(:,i,1).*CDWGB(:,j,1)...
                      + PermK(:,1,2).*CDWGB(:,i,1).*CDWGB(:,j,2)...
                      + PermK(:,2,1).*CDWGB(:,i,2).*CDWGB(:,j,1)...
                      + PermK(:,2,2).*CDWGB(:,i,2).*CDWGB(:,j,2)).*area;
  end
end

%% Computing EltStabMat: Elementwise Stabilizer Matrix 
EltStabMat = zeros(NumEms,7,7);
% Term #1 for edge #1 
ESM = zeros(NumEms,7,7);
leneg = sqrt((x1-x2).^2+(y1-y2).^2);
dg = dm1.*leneg;
Xm = 0.5*(x1+x2) - xc;
Ym = 0.5*(y1+y2) - yc;
ESM(:,1,1) =  dg;  ESM(:,2,2) = dg.*(Xm.^2);  ESM(:,3,3) = dg.*(Ym.^2);
ESM(:,1,2) =  dg.*Xm;      ESM(:,2,1) = ESM(:,1,2);
ESM(:,1,3) =  dg.*Ym;      ESM(:,3,1) = ESM(:,1,3);
ESM(:,2,3) =  dg.*Xm.*Ym;  ESM(:,3,2) = ESM(:,2,3);
ESM(:,1,4) = -dg;          ESM(:,4,1) = ESM(:,1,4);
ESM(:,2,4) = -dg.*Xm;      ESM(:,4,2) = ESM(:,2,4);
ESM(:,3,4) = -dg.*Ym;      ESM(:,4,3) = ESM(:,3,4);
ESM(:,4,4) =  dg;
EltStabMat = EltStabMat + ESM;
% Term #2 for edge #2 
ESM = zeros(NumEms,7,7);
leneg = sqrt((x2-x3).^2+(y2-y3).^2);
dg = dm1.*leneg;
Xm = 0.5*(x2+x3) - xc;
Ym = 0.5*(y2+y3) - yc;
ESM(:,1,1) =  dg;  ESM(:,2,2) = dg.*(Xm.^2);  ESM(:,3,3) = dg.*(Ym.^2);
ESM(:,1,2) =  dg.*Xm;      ESM(:,2,1) = ESM(:,1,2);
ESM(:,1,3) =  dg.*Ym;      ESM(:,3,1) = ESM(:,1,3);
ESM(:,2,3) =  dg.*Xm.*Ym;  ESM(:,3,2) = ESM(:,2,3);
ESM(:,1,5) = -dg;          ESM(:,5,1) = ESM(:,1,5);
ESM(:,2,5) = -dg.*Xm;      ESM(:,5,2) = ESM(:,2,5);
ESM(:,3,5) = -dg.*Ym;      ESM(:,5,3) = ESM(:,3,5);
ESM(:,5,5) =  dg;
EltStabMat = EltStabMat + ESM;
% Term #3 for edge #3 
ESM = zeros(NumEms,7,7);
leneg = sqrt((x3-x4).^2+(y3-y4).^2);
dg = dm1.*leneg;
Xm = 0.5*(x3+x4) - xc;
Ym = 0.5*(y3+y4) - yc;
ESM(:,1,1) =  dg;  ESM(:,2,2) = dg.*(Xm.^2);  ESM(:,3,3) = dg.*(Ym.^2);
ESM(:,1,2) =  dg.*Xm;      ESM(:,2,1) = ESM(:,1,2);
ESM(:,1,3) =  dg.*Ym;      ESM(:,3,1) = ESM(:,1,3);
ESM(:,2,3) =  dg.*Xm.*Ym;  ESM(:,3,2) = ESM(:,2,3);
ESM(:,1,6) = -dg;          ESM(:,6,1) = ESM(:,1,6);
ESM(:,2,6) = -dg.*Xm;      ESM(:,6,2) = ESM(:,2,6);
ESM(:,3,6) = -dg.*Ym;      ESM(:,6,3) = ESM(:,3,6);
ESM(:,6,6) =  dg;
EltStabMat = EltStabMat + ESM;
% Term #4 for edge #4 
ESM = zeros(NumEms,7,7);
leneg = sqrt((x4-x1).^2+(y4-y1).^2);
dg = dm1.*leneg;
Xm = 0.5*(x4+x1) - xc;
Ym = 0.5*(y4+y1) - yc;
ESM(:,1,1) =  dg;  ESM(:,2,2) = dg.*(Xm.^2);  ESM(:,3,3) = dg.*(Ym.^2);
ESM(:,1,2) =  dg.*Xm;      ESM(:,2,1) = ESM(:,1,2);
ESM(:,1,3) =  dg.*Ym;      ESM(:,3,1) = ESM(:,1,3);
ESM(:,2,3) =  dg.*Xm.*Ym;  ESM(:,3,2) = ESM(:,2,3);
ESM(:,1,7) = -dg;          ESM(:,7,1) = ESM(:,1,7);
ESM(:,2,7) = -dg.*Xm;      ESM(:,7,2) = ESM(:,2,7);
ESM(:,3,7) = -dg.*Ym;      ESM(:,7,3) = ESM(:,3,7);
ESM(:,7,7) =  dg;
EltStabMat = EltStabMat + ESM;

%% NOTE: FOR CONVENIENCE, BUT USE TOO MUCH MEMORY 
%% Elementwise matrices 
EltMat = EltStfMat + rho*EltStabMat;

%% Assembling GlbMat 
% Element interior interaction 
for i=1:3
  II = 3*(0:NumEms-1) + i;
  for j=1:3
    JJ = 3*(0:NumEms-1) + j;
    GlbMat = GlbMat + sparse(II,JJ,EltMat(:,i,j),DOFs,DOFs);    
  end
end
% Element-edge interaction 
for i=1:3
  II = 3*(0:NumEms-1) + i;
  for j=1:4
    JJ = 3*NumEms + QuadriMesh.elem2edge(:,j);
    GlbMat = GlbMat + sparse(II,JJ,EltMat(:,i,3+j),DOFs,DOFs);
    GlbMat = GlbMat + sparse(JJ,II,EltMat(:,i,3+j),DOFs,DOFs);    
  end
end
% Edge-edge interaction 
for i=1:4
  II = 3*NumEms + QuadriMesh.elem2edge(:,i);
  for j=i:4  % Utilizing symmetry 
    JJ = 3*NumEms + QuadriMesh.elem2edge(:,j);
    GlbMat = GlbMat + sparse(II,JJ,EltMat(:,3+i,3+j),DOFs,DOFs);
    if (j>i) 
      GlbMat = GlbMat + sparse(JJ,II,EltMat(:,3+i,3+j),DOFs,DOFs);
    end
  end
end

%% Assembling GlbRHS 
GlbRHS = zeros(DOFs,1);
tmp1 = zeros(NumEms,1);  tmpX = zeros(NumEms,1);  tmpY = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,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);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  fval = EqnBC.fxnf(qp);
  tmp1  = tmp1  + GAUSSQUAD.RECT(k,3)*jac.*fval;
  tmpX  = tmpX  + GAUSSQUAD.RECT(k,3)*jac.*(fval.*X);
  tmpY  = tmpY  + GAUSSQUAD.RECT(k,3)*jac.*(fval.*Y);
end
GlbRHS(3*(0:NumEms-1)+1) = GlbRHS(3*(0:NumEms-1)+1) + tmp1;
GlbRHS(3*(0:NumEms-1)+2) = GlbRHS(3*(0:NumEms-1)+2) + tmpX;
GlbRHS(3*(0:NumEms-1)+3) = GlbRHS(3*(0:NumEms-1)+3) + tmpY;

%% Incorporating boundary conditions: Neumann as natural 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  x1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdge,2),2);
  LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  qp = zeros(NumNeumannEgs,2);
  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;
    GlbRHS(NumEms+NeumannEdge) = GlbRHS(NumEms+NeumannEdge)...
      - GAUSSQUAD.LINE(k,3) * EqnBC.fxnuN(qp) .* LenNeumannEg;
  end
end

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

%% Incorporating boundary conditions: Dirichlet as essential 
% Approach I: working on all these edges simultaneously 
if NumDirichletEgs>0 
  x1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdge,2),2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  qp = zeros(NumDirichletEgs,2);
  pDavg = zeros(NumDirichletEgs,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;
    pDavg = pDavg + GAUSSQUAD.LINE(k,3)*EqnBC.fxnpD(qp);
  end
  sln(3*NumEms+DirichletEdge) = pDavg;
end

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

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

return;