function sln = SteadyCD_WG_RectQ0Q0RT0_AsmSlv(...
  EqnBC,RectMesh,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Steady conv.-diff.eqn.: WG(Q0,Q0;RT[0]) on a rect.mesh: Assembly & Solving
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;
% Assuming RectMesh.flag>=2 
area = RectMesh.area;
e2g = RectMesh.elem2edge;

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

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

%% For convection: 
% Note the negative sign !!
% ArrayConv = zeros(NumEms,5);
ArrayConv = WG_RectQ0Q0RT0_MixProdVecGradValBasFxn(EqnBC.fxnv, ...
  RectMesh,GAUSSQUAD);
II = (1:NumEms);  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,-ArrayConv(:,1),DOFs,DOFs);
JJ = (1:NumEms);
for i=1:4
  II = NumEms + e2g(:,i);
  GlbMat = GlbMat + sparse(II,JJ,-ArrayConv(:,i+1),DOFs,DOFs);
end

%% For diffusion: 
D = EqnBC.D;
DiffD = zeros(NumEms,2,2);
DiffD(:,1,1) = D;
DiffD(:,2,2) = D;
GMK = Hdiv_RectRT0_NmlzBas_GramMatK(RectMesh,DiffD);
CDWGB = WG_RectQ0Q0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(RectMesh);
ArrayEE = zeros(NumEms,1);
ArrayEG = zeros(NumEms,4);
ArrayGG = zeros(NumEms,4,4);
%% JL20151225: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEms
  EltGM = squeeze(GMK(ie,:,:));
  ArrayEE(ie) = squeeze(CDWGB(ie,1,:))' * EltGM * squeeze(CDWGB(ie,1,:));
  for i=1:4
    ArrayEG(ie,i) = squeeze(CDWGB(ie,1,:))' * EltGM * squeeze(CDWGB(ie,1+i,:));
    for j=1:4
      ArrayGG(ie,i,j) = squeeze(CDWGB(ie,1+i,:))' * EltGM * squeeze(CDWGB(ie,1+j,:));
    end
  end
end
% Part 1: element-element interaction 
II = (1:NumEms);  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,ArrayEE,DOFs,DOFs);
% Part 2: element-edge interaction 
II = (1:NumEms);
for j=1:4
  JJ = NumEms + e2g(:,j);
  GlbMat = GlbMat + sparse(II,JJ,ArrayEG(:,j),DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,ArrayEG(:,j),DOFs,DOFs);
end
% Part 3: edge-edge interaction 
for i=1:4
  II = NumEms + e2g(:,i);
  for j=i:4  % Utilizing symmetry 
    JJ = NumEms + e2g(:,j);
    if (j==i) 
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,i),DOFs,DOFs);
    else
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
      GlbMat = GlbMat + sparse(JJ,II,ArrayGG(:,i,j),DOFs,DOFs);
    end
  end
end

%% Assembling GlbRHS 
% Applying a Gaussian quadrature on rectangular elements 
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);
GlbRHS = zeros(DOFs,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.RECT(k,1)*x1 + GAUSSQUADRECT1(k,1)*x2;
  qp(:,2) = GAUSSQUAD.RECT(k,2)*y1 + GAUSSQUADRECT1(k,2)*y2;
  GlbRHS(1:NumEms) = GlbRHS(1:NumEms) + GAUSSQUAD.RECT(k,3)*EqnBC.fxnf(qp);
end
GlbRHS(1:NumEms) = GlbRHS(1:NumEms).*area;

%% Incorporating boundary conditions: Neumann as natural 
% Approach I: working on all these edges simultaneously 
if NumNeumannEgs>0
  x1 = RectMesh.node(RectMesh.edge(NeumannEdge,1),1);
  y1 = RectMesh.node(RectMesh.edge(NeumannEdge,1),2);
  x2 = RectMesh.node(RectMesh.edge(NeumannEdge,2),1);
  y2 = RectMesh.node(RectMesh.edge(NeumannEdge,2),2);
  LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
  % LenNeumannEg = norm(RectMesh.node(RectMesh.edge(NeumannEdge,2),:)...
  %   - RectMesh.node(RectMesh.edge(NeumannEdge,1),:),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.fxncN(qp).*LenNeumannEg;
  end
end

%% For reducing... 
% disp('Adjusting the global linear system...'); 
flag = zeros(DOFs,1);
flag(NumEms+DirichletEdge) = ones(NumDirichletEgs,1);
EmFreeEg = find(~flag);
% sln vector to be used 
sln = zeros(DOFs,1);

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

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

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

return;