function sln = TransCDRnv_IE_WG_RectQ0Q0RT0_AsmSlv(EqnBIC,CofRT0NB,...
  RectMesh,tyme,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Solving transient convection-diffusion-reaction eqn. (numer. velocity)
% by Implicit Euler + WG(Q0,Q0;RT[0]) on a rectangular mesh 
% Assuming a constant reaction coefficient
% James Liu, ColoState; 2012/07--2017/02

%% EqnBIC 
D = EqnBIC.D;

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

%% Tyme info
NT = tyme.NT;  MT = tyme.MT;  Deltat = tyme.Deltat;

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

%% Setup 
DOFs = NumEms + NumEgs;
% GlbMat = sparse(DOFs,DOFs);

%% Assembling global matrices: 
% Four parts: mass, convection, diffusion, reaction

%% Assembling global matrices: For mass:
GlbMatMass = sparse(DOFs,DOFs);
II = (1:NumEms);  JJ = II;
GlbMatMass = GlbMatMass + sparse(II, JJ, RectMesh.area, DOFs, DOFs);

%% Assembling global matrices: For convection:
GlbMatConv = sparse(DOFs,DOFs);
% Assuming numerical velocity in RT[0] 
GMRT0 = Hdiv_RectRT0_NmlzBas_GramMat(RectMesh);
CDWGB = WG_RectQ0Q0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(RectMesh);
ArrayConv = zeros(NumEms,5);
for i=1:5
  ArrayConv(:,i) = sum(CofRT0NB.*squeeze(CDWGB(:,i,:)).*GMRT0, 2);
end
% Assuming numerical velocity in RT[1]: To be added later  
% Now assembly: Note the negative sign !!
II = (1:NumEms);  JJ = II;
GlbMatConv = GlbMatConv + sparse(II,JJ,-ArrayConv(:,1),DOFs,DOFs);
JJ = (1:NumEms);
for i=1:4
  II = NumEms + RectMesh.elem2edge(:,i);
  GlbMatConv = GlbMatConv + sparse(II,JJ,-ArrayConv(:,i+1),DOFs,DOFs);
end

%% Assembling global matrices: For diffusion:
GlbMatDiff = sparse(DOFs,DOFs);
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;
GlbMatDiff = GlbMatDiff + sparse(II,JJ,ArrayEE,DOFs,DOFs);
% Part 2: element-edge interaction 
II = (1:NumEms);
for j=1:4
  JJ = NumEms + RectMesh.elem2edge(:,j);
  GlbMatDiff = GlbMatDiff + sparse(II,JJ,ArrayEG(:,j),DOFs,DOFs);
  GlbMatDiff = GlbMatDiff + sparse(JJ,II,ArrayEG(:,j),DOFs,DOFs);
end
% Part 3: edge-edge interaction 
for i=1:4
  II = NumEms + RectMesh.elem2edge(:,i);
  for j=i:4  % Utilizing symmetry 
    JJ = NumEms + RectMesh.elem2edge(:,j);
    if (j==i) 
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,ArrayGG(:,i,i),DOFs,DOFs);
    else
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
      GlbMatDiff = GlbMatDiff + sparse(JJ,II,ArrayGG(:,i,j),DOFs,DOFs);
    end
  end
end

%% Assembling global matrices: For (linear) reaction:
GlbMatReac = sparse(DOFs,DOFs);
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);
ArrayReacEE = zeros(NumEms,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;
  ArrayReacEE = ArrayReacEE + GAUSSQUAD.RECT(k,3)*EqnBIC.fxnR(qp);
end
ArrayReacEE = ArrayReacEE .* RectMesh.area;
II = (1:NumEms);  JJ = II;
GlbMatReac = GlbMatReac + sparse(II, JJ, ArrayReacEE, DOFs, DOFs);

%% Assembling global matrices: Combined
GlbMat = GlbMatMass + Deltat * (GlbMatConv + GlbMatDiff + GlbMatReac);
clear  GlbMatConv  GlbMatDiff  GlbMatReac;

%% Initial numerical solution 
sold = zeros(DOFs,1);
% For element interior 
sold(1:NumEms) = zeros(NumEms,1);
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);
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;
  sold(1:NumEms) = sold(1:NumEms) ...
    + GAUSSQUAD.RECT(k,3)*EqnBIC.fxnc0(qp);  % */area 
end
% For mesh skeleton 
x1 = RectMesh.node(RectMesh.edge(:,1),1);
y1 = RectMesh.node(RectMesh.edge(:,1),2);
x2 = RectMesh.node(RectMesh.edge(:,2),1);
y2 = RectMesh.node(RectMesh.edge(:,2),2);
NumQuadPts = size(GAUSSQUAD.LINE,1);
qp = zeros(NumEgs,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;
  sold(NumEms+(1:NumEgs)) = sold(NumEms+(1:NumEgs)) ... 
    + GAUSSQUAD.LINE(k,3)*EqnBIC.fxnc0(qp);  % */ LenEg
end 
% Storing 
sln(:,1) = sold;

%% Time-marching 
for n=1:MT
  % tnew = n * Deltat;
  tnew = tyme.tm(n+1);

  % For the global right-hand side 
  GlbRHS = zeros(DOFs,1);

  % Assembling the source term at the new time 
  EltSrc = zeros(NumEms,1);  
  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);
  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;
    EltSrc = EltSrc + GAUSSQUAD.RECT(k,3) * EqnBIC.fxnf(qp,tnew);
  end
  EltSrc = Deltat * (EltSrc.*RectMesh.area) ;
  GlbRHS(1:NumEms,1) = GlbRHS(1:NumEms,1) + EltSrc;

  % Assembling the old mass 
  GlbRHS = GlbRHS + GlbMatMass*sold;

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

  % Incorporating boundary conditions: Dirichlet as essential 
  % Approach I: working on all these edges simultaneously 
  if NumDirichletEgs>0 
    cDavg = zeros(NumDirichletEgs,1);
    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);
    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)*EqnBIC.fxncD(qp,tnew);
    end
    snew(NumEms+DirichletEdge) = cDavg;
  end

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

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

  % Storing  
  sln(:,n+1) = snew;

  % Marching to the next time step  
  sold = snew;
end

return;