function sln = Parab_CN_CG_RectQ1_AsmSlv(...
  EqnBIC,RectMesh,tyme,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Solving the parabolic equation 
% by Crank-Nicolson + Continuous Galerkin Q1 on a rectangular mesh
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info
NumNds = RectMesh.NumNds;
NumEms = RectMesh.NumEms;
area = RectMesh.area;

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

%% Boundary info 
NumDirichletEgs = size(DirichletEdge,1);
NumNeumannEgs = size(NeumannEdge,1);
% 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 "static" global matrix: Mass 
% For mass: elementwise 
% EltMatMass = zeros(NumEms,3,3);
EltMatMass = CG_RectQ1_GramMat(RectMesh);
% For mass: global 
GlbMatMass = sparse(DOFs,DOFs);
for i=1:4
  II = RectMesh.elem(:,i);
  for j=1:4
    JJ = RectMesh.elem(:,j);
    GlbMatMass = GlbMatMass + sparse(II,JJ,EltMatMass(:,i,j),DOFs,DOFs);
  end
end

%% Assembling the "static" global matrix: Stiffness/Diffusion  
% For diffusion: elementwise 
% EltMatDiff = zeros(NumEms,4,4);
EltMatDiff = CG_RectQ1_GradGradMat(RectMesh); 
% For diffusion: global 
GlbMatDiff = sparse(DOFs,DOFs);
for i=1:4
  II = RectMesh.elem(:,i);
  for j=i:4  % Utilizing symmetry
    JJ = RectMesh.elem(:,j);
    if (j==i) 
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,EltMatDiff(:,i,j),DOFs,DOFs);
    else
      GlbMatDiff = GlbMatDiff + sparse(II,JJ,EltMatDiff(:,i,j),DOFs,DOFs);
      GlbMatDiff = GlbMatDiff + sparse(JJ,II,EltMatDiff(:,i,j),DOFs,DOFs);
    end
  end
end

%% Assembling the "static" global matrix: Combined 
GlbMatL = GlbMatMass + 0.5*Deltat * GlbMatDiff;  % matrix on the left-side 
GlbMatR = GlbMatMass - 0.5*Deltat * GlbMatDiff;  % matrix on the right-side
clear  GlbMatDiff;

%% Initial numerical solution
sold = zeros(DOFs,1);
% Setting the initial numerical solution by interpolation 
sold = EqnBIC.fxnc0(RectMesh.node);
% Storing 
sln(:,1) = sold;

% For convenience
k1 = RectMesh.elem(:,1);  k2 = RectMesh.elem(:,3);
x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2);
x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);

% Assembling the source term at the initial time 
NdSrcOld = zeros(NumNds,1);
EltSrc = zeros(NumEms,4);
qp = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
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 = EqnBIC.fxnf(qp,tyme.tm(1));
  EltSrc(:,1) = EltSrc(:,1) + GAUSSQUAD.RECT(k,3) * fval .* (1-X).*(1-Y);
  EltSrc(:,2) = EltSrc(:,2) + GAUSSQUAD.RECT(k,3) * fval .*     X.*(1-Y);
  EltSrc(:,3) = EltSrc(:,3) + GAUSSQUAD.RECT(k,3) * fval .*     X.*Y;
  EltSrc(:,4) = EltSrc(:,4) + GAUSSQUAD.RECT(k,3) * fval .* (1-X).*Y;
end
EltSrc = EltSrc .* [area,area,area,area];
% JL20160321: TO BE REVISED FOR EFFECIENCY 
% Converting 
for ie=1:NumEms
  for j=1:4
    NdSrcOld(RectMesh.elem(ie,j)) = ...
      NdSrcOld(RectMesh.elem(ie,j)) + EltSrc(ie,j);
  end
end

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

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

  % Assembling the source term at the new time 
  NdSrcNew = zeros(NumNds,1);
  EltSrc = zeros(NumEms,4);
  qp = zeros(NumEms,2);
  NumQuadPts = size(GAUSSQUAD.RECT,1);
  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 = EqnBIC.fxnf(qp,tnew);
    EltSrc(:,1) = EltSrc(:,1) + GAUSSQUAD.RECT(k,3) * fval .* (1-X).*(1-Y);
    EltSrc(:,2) = EltSrc(:,2) + GAUSSQUAD.RECT(k,3) * fval .*     X.*(1-Y);
    EltSrc(:,3) = EltSrc(:,3) + GAUSSQUAD.RECT(k,3) * fval .*     X.*Y;
    EltSrc(:,4) = EltSrc(:,4) + GAUSSQUAD.RECT(k,3) * fval .* (1-X).*Y;
  end
  EltSrc = EltSrc .* [area,area,area,area];
  % JL20160321: TO BE REVISED FOR EFFECIENCY 
  % Converting 
  for ie=1:NumEms
    for j=1:4
      NdSrcNew(RectMesh.elem(ie,j)) = ...
        NdSrcNew(RectMesh.elem(ie,j)) + EltSrc(ie,j);
    end
  end

  % Assembling the old and new source terms
  GlbRHS = Deltat * 0.5*(NdSrcOld+NdSrcNew);

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

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

  % Incorporating boundary conditions: Dirichlet as essential 
  % Approach I: working on all these nodes simultaneously 
  % Assuming NumDirichletNds>0
  DiriNd = RectMesh.node(DirichletNode,:);
  cDval = EqnBIC.fxncD(DiriNd,tnew);
  snew(DirichletNode) = cDval;

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

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

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

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

return;