%% JL20170206: Old code for 2-phase problems, TO BE REVISED !
function [sln,GlbMat,GlbRHS] = Darcy_WG_RectQ0P0RT0s_AsmSlv(EqnBC,...
  RectMesh,PermKs,DirichletEdge,NeumannEdge,GAUSSQUAD)
%% Darcy: Weak Galerkin (Q0,Q0;RT[0]) on rectangular meshes 
% Assuming PermK is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2017/02

%% Mesh info 
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;

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

%% Auxiliary arrays 
GMRT0 = Hdiv_RectRT0_NmlzBas_GramMat(RectMesh);
CDWGB = WG_RectQ0Q0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(RectMesh);

%% Setting up 
DOFs = NumEms + NumEgs;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);
% sln = zeros(DOFs,1);

%% Assembling the global coefficient matrix GlbMat 
GlbMat = sparse(DOFs,DOFs);
% Part I: element-element interaction 
Aij = PermKs .* sum(squeeze(CDWGB(:,1,:)).*GMRT0.*squeeze(CDWGB(:,1,:)),2);
II = (1:NumEms);  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,Aij,DOFs,DOFs);
% Part II: element-edge interaction 
for j=1:4
  Aij = PermKs .* sum(squeeze(CDWGB(:,1,:)).*GMRT0.*squeeze(CDWGB(:,j+1,:)),2);
  II = (1:NumEms);  JJ = NumEms + RectMesh.elem2edge(:,j);
  MatB = sparse(II,JJ,Aij,DOFs,DOFs);
  GlbMat = GlbMat + MatB + MatB';
end
% Part III: edge-edge interaction 
for i=1:4
  II = NumEms + RectMesh.elem2edge(:,i);
  for j=i:4  % Utilizing symmetry 
    JJ = NumEms + RectMesh.elem2edge(:,j);
    Aij = PermKs .* sum(squeeze(CDWGB(:,i+1,:)).*GMRT0.*squeeze(CDWGB(:,j+1,:)),2);
    if (j==i) 
      GlbMat = GlbMat + sparse(II,JJ,Aij,DOFs,DOFs);
    else
      GlbMat = GlbMat + sparse([II,JJ],[JJ,II],[Aij;Aij],DOFs,DOFs);
    end
  end
end
clear  PermK  Aij  MatB;

%% 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) .* RectMesh.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.fxnuN(qp) .* LenNeumannEg;
  end
end
% % Approach II: (the old-style) loop thru these edges 
% if NumNeumannEgs>0
%   for ig=1:NumNeumannEgs
%     LblEg = NeumannEdge(ig);
%     LblVrtx = RectMesh.edge(LblEg,1:2);
%     vrtx = RectMesh.node(LblVrtx,1:2);
%     LenEg = norm(vrtx(2,:)-vrtx(1,:),2);
%     QuadPt = GAUSSQUAD.LINE(:,1:2)*vrtx;
%     uNval = EqnBC.fxnuN(QuadPt);
%     pos = NumEms + LblEg;
%     GlbRHS(pos) = GlbRHS(pos) - dot(GAUSSQUAD.LINE(:,3),uNval)*LenEg;
%   end
% end

%% For reducing... 
% disp('Adjusting the global linear system...'); 
flag = zeros(DOFs,1);
flag(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 = 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);
  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(NumEms+DirichletEdge) = pDavg;
end
% % Approach II: (the old-style) loop thru these edges 
% if NumDirichletEgs>0 
%   for ig=1:NumDirichletEgs
%     LblEg = DirichletEdge(ig);
%     LblVrtx = RectMesh.edge(LblEg,1:2);
%     vrtx = RectMesh.node(LblVrtx,1:2);
%     QuadPt = GAUSSQUAD.LINE(:,1:2)*vrtx;
%     pDval = EqnBC.fxnpD(QuadPt);
%     pDavg = dot(GAUSSQUAD.LINE(:,3),pDval);  % *LenEg/LenEg
%     sln(NumEms+LblEg) = pDavg;
%   end
% end

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

% %% Solving the full global linear system directly 
% disp('Solving the full global linear system directly...'); 
% % First modifying the global linear system 
% RHS = GlbRHS;
% GlbRHS(NumEms+DirichletEdge) = sln(NumEms+DirichletEdge);
% GM = sparse(DOFs,DOFs);
% GM(EmFreeEg,EmFreeEg) = GlbMat(EmFreeEg,EmFreeEg);
% GM(NumEms+DirichletEdge,NumEms+DirichletEdge) = eye(NumDirichletEgs);
% % Then solving 
% sln  = GM\RHS;

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

%% JL20130728: BEING ADDED: SOLVING BY SCHUR COMPLEMENT 
% Note: No interaction across elements, Amm and its inverse are digonal 
% MGM = GlbMat(EmFreeEg,EmFreeEg);  % Squeezed (Modified Global Matrix)
% RHS = GlbRHS(EmFreeEg);  % Squeezed (Modified RHS)
% A11i = spdiags(1./diag(MGM(1:NumEms,1:NumEms)),0,NumEms,NumEms);
% A12 = MGM(1:NumEms,NumEms+1:end); 
% A21 = A12';
% A22 = MGM(NumEms+1:end,NumEms+1:end);
% A = A22 - A21*A11i*A12;
% b1 = RHS(1:NumEms);
% b2 = RHS(NumEms+1,end);
% b = b2 - A21*A11i*b1;
% % Solving lin. sys.: A z = b 
% sln(1:NumEms) = A11i*(b1-A12*z); 
% sln(NumEms+1:end) = z "expanded" 

return;