function [sln,GlbMat,GlbRHS] = Darcy_CG_RectQ1s_AsmSlv(...
  EqnBC,RectMesh,PermKs,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Darcy: CG: Q1 on a rectangular mesh: Assembling & Solving 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2016/06 

%% Mesh info 
NumNds = RectMesh.NumNds;
NumEms = RectMesh.NumEms;
% Assuming TriMesh.flag>=2 
% area = RectMesh.area;

%% Boundary info 
NumDirichletEgs = size(DirichletEdge,1);
NumNeumannEgs = size(NeumannEdge,1);
% JL20151229: TO BE REVISED FOR EFFICIENCY 
DirichletNodeFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgs
  k1 = RectMesh.edge(DirichletEdge(ig),1);
  k2 = RectMesh.edge(DirichletEdge(ig),2);
  DirichletNodeFlag(k1) = 1;
  DirichletNodeFlag(k2) = 1;
end
DirichletNode = find(DirichletNodeFlag);
NumDirichletNds = sum(DirichletNodeFlag);

%% Setting up the global linear system 
DOFs = NumNds;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);
% sln = zeros(DOFs,1);

%% Assembling the global coefficient matrix GlbMat 
% Assuming PermK is an elementwise scalar 
EGGM = CG_RectQ1_GradGradMat(RectMesh);
GlbMat = sparse(DOFs,DOFs);
for i=1:4
  II = RectMesh.elem(:,i);
  for j=1:4
    JJ = RectMesh.elem(:,j);
    Aij = PermKs .* EGGM(:,i,j);
    GlbMat = GlbMat + sparse(II,JJ,Aij,NumNds,NumNds);
  end
end

%% JL20131004: Old style temporarily 
%% Assembling GlbRHS: The source term 
GlbRHS = zeros(DOFs,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
% tmp = zeros(4,1);
qp = zeros(1,2);
for ie=1:NumEms
  k1 = RectMesh.elem(ie,1);  k2 = RectMesh.elem(ie,3);
  x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2); 
  x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);
  ar = (x2-x1) * (y2-y1);
  % tmp = GlbRHS(RectMesh.elem(ie,1:4));
  for k=1:NumQuadPts
    qp(1) = GAUSSQUAD.RECT(k,1)*x1 + (1-GAUSSQUAD.RECT(k,1))*x2;
    qp(2) = GAUSSQUAD.RECT(k,2)*y1 + (1-GAUSSQUAD.RECT(k,2))*y2;
    X = (qp(1)-x1)/(x2-x1);  Y = (qp(2)-y1)/(y2-y1);    
    fval = EqnBC.fxnf(qp);
    GlbRHS(RectMesh.elem(ie,1)) = GlbRHS(RectMesh.elem(ie,1)) ...
      + GAUSSQUAD.RECT(k,3) * ar * fval * (1-X)*(1-Y);
    GlbRHS(RectMesh.elem(ie,2)) = GlbRHS(RectMesh.elem(ie,2)) ...
      + GAUSSQUAD.RECT(k,3) * ar * fval *     X*(1-Y);
    GlbRHS(RectMesh.elem(ie,3)) = GlbRHS(RectMesh.elem(ie,3)) ...
      + GAUSSQUAD.RECT(k,3) * ar * fval *     X*Y;
    GlbRHS(RectMesh.elem(ie,4)) = GlbRHS(RectMesh.elem(ie,4)) ...
      + GAUSSQUAD.RECT(k,3) * ar * fval * (1-X)*Y;
  end
end

%% Incorporating boundary conditions: Neumann as natural 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  k1 = RectMesh.edge(NeumannEdge,1);  
  k2 = RectMesh.edge(NeumannEdge,2);
  x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2);
  x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);
  % LenEgNeum = 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;
    phi1 = sqrt((x2-qp(:,1)).^2+(y2-qp(:,2)).^2);  % ./LenEgNeum;
    phi2 = sqrt((qp(:,1)-x1).^2+(qp(:,2)-y1).^2);  % ./LenEgNeum;
    uNval = EqnBC.fxnuN(qp);
    GlbRHS(k1) = GlbRHS(k1) - GAUSSQUAD.LINE(k,3)*(uNval.*phi1);
    GlbRHS(k2) = GlbRHS(k2) - GAUSSQUAD.LINE(k,3)*(uNval.*phi2);
  end
end

%% For reducing... 
% disp('Adjusting the global linear system...'); 
flag = zeros(DOFs,1);
flag(DirichletNode) = ones(NumDirichletNds,1);
FreeNd = find(~flag);
% sln vector to be used 
sln = zeros(DOFs,1);

%% Incorporating boundary conditions: Dirichlet as essential 
% Approach I: working on all these nodes simultaneously 
% Assuming NumDirichletNds>0
DiriNd = RectMesh.node(DirichletNode,:);
pDval = EqnBC.fxnpD(DiriNd);
sln(DirichletNode) = pDval;

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

%% Solving the reduced global linear system directly 
% disp('Solving the reduced global linear system directly...'); 
sln(FreeNd)  = GlbMat(FreeNd,FreeNd) \ GlbRHS(FreeNd);

return;