function [CofNS,GlbMat,GlbRHS] = Darcy_DG_TriP1_AsmSlv(EqnBC,...
  TriMesh,PermK,DirichletEdge,NeumannEdge,GAUSSQUAD,alpha,beta)
%% Darcy: DG: TriP1: Assembling & solving 
% Assuming PermK is an elementwise scalar constant 
% James Liu, ColoState; 2012/07--2016/06 

%% Mesh info 
NumEms = TriMesh.NumEms;
% NumEgs = TriMesh.NumEgs;
area = TriMesh.area;
LenEg = TriMesh.LenEg;

%% Boundary info 
NumDiriEgs = size(DirichletEdge,1);
NumNeumEgs = size(NeumannEdge,1);
NumIntrEgs = sum(TriMesh.BndryEdge==0);
IntrEdge = find(TriMesh.BndryEdge==0);

% %% Penalty factor and form parameter 
% % a magic number as suggested by Ohannes Karakashian 
% alpha = 6;
% Form parameter !!!
% beta = 1;  % for nonsymmetric formulation 

%% Setting up the global linear system 
DOFs = 3*NumEms;
GlbMat = sparse(DOFs,DOFs);
GlbRHS = zeros(NumEms,3);  % This is nonconventional! 

%% Assembling the global coeff. matrix GlbMat 

%% GlbMat: Term 1: grad-grad on all triangular elements 
MatGG = zeros(NumEms,3,3);
MatGG(:,2,2) = area;  
MatGG(:,3,3) = area;
% Assuming PermK is an elementwise scalar constant 
for i=2:3
  II = 3*(0:NumEms-1) + i;
  Aii = PermK.*MatGG(:,i,i);
  GlbMat = GlbMat + sparse(II,II,Aii,DOFs,DOFs);
end
clear MatGG Aii;

%% JL20131016: TO BE REVISED FOR EFFECIENCY 
%% GlbMat: Term 2&3: Terms from integration by parts for interior edges 
A11 = zeros(NumIntrEgs,3,3);  A12 = zeros(NumIntrEgs,3,3);  
A21 = zeros(NumIntrEgs,3,3);  A22 = zeros(NumIntrEgs,3,3);  
B11 = zeros(NumIntrEgs,3,3);  B12 = zeros(NumIntrEgs,3,3);  
B21 = zeros(NumIntrEgs,3,3);  B22 = zeros(NumIntrEgs,3,3);  
GradBasFxn = [0,0; 1,0; 0,1];
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig=1:NumIntrEgs
  x1 = TriMesh.node(TriMesh.edge(IntrEdge(ig),1),1);
  y1 = TriMesh.node(TriMesh.edge(IntrEdge(ig),1),2);
  x2 = TriMesh.node(TriMesh.edge(IntrEdge(ig),2),1);
  y2 = TriMesh.node(TriMesh.edge(IntrEdge(ig),2),2);
  nml = repmat(TriMesh.NmlEg(IntrEdge(ig),:),3,1);
  ie1 = TriMesh.edge2elem(IntrEdge(ig),1);
  ie2 = TriMesh.edge2elem(IntrEdge(ig),2);
  PermK1 = PermK(ie1);  PermK2 = PermK(ie2);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    X = qp(1);  Y = qp(2);
    BasFxn = [1; X; Y];
    Kdpgn1 = PermK1*dot(GradBasFxn,nml,2);
    Kdpgn2 = PermK2*dot(GradBasFxn,nml,2);
    for i=1:3
      for j=1:3
        A11(ig,i,j) = A11(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn(i)*Kdpgn1(j)*LenEg(IntrEdge(ig));
        A12(ig,i,j) = A12(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn(i)*Kdpgn2(j)*LenEg(IntrEdge(ig));
        A21(ig,i,j) = A21(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn(i)*Kdpgn1(j)*LenEg(IntrEdge(ig));
        A22(ig,i,j) = A22(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn(i)*Kdpgn2(j)*LenEg(IntrEdge(ig));
      end
    end
  end
  % JL20131016: NEEDS DOUBLE-CHECK !!!
  A11(ig,:,:) =  0.5*A11(ig,:,:);
  A12(ig,:,:) =  0.5*A12(ig,:,:);
  A21(ig,:,:) = -0.5*A21(ig,:,:);
  A22(ig,:,:) = -0.5*A22(ig,:,:);
  B11(ig,:,:) = squeeze(A11(ig,:,:))';  
  B12(ig,:,:) = squeeze(A21(ig,:,:))';
  B21(ig,:,:) = squeeze(A12(ig,:,:))';
  B22(ig,:,:) = squeeze(A22(ig,:,:))';
end
% Now assembling 
LblElt1 = TriMesh.edge2elem(IntrEdge,1);
LblElt2 = TriMesh.edge2elem(IntrEdge,2);
for i=1:3
  II1 = 3*(LblElt1-1) + i;
  II2 = 3*(LblElt2-1) + i;
  for j=1:3
    JJ1 = 3*(LblElt1-1) + j;
    JJ2 = 3*(LblElt2-1) + j;
    GlbMat = GlbMat + sparse(II1,JJ1,-A11(:,i,j)+beta*B11(:,i,j),...
      DOFs,DOFs);
    GlbMat = GlbMat + sparse(II1,JJ2,-A12(:,i,j)+beta*B12(:,i,j),...
      DOFs,DOFs);
    GlbMat = GlbMat + sparse(II2,JJ1,-A21(:,i,j)+beta*B21(:,i,j),...
      DOFs,DOFs);
    GlbMat = GlbMat + sparse(II2,JJ2,-A22(:,i,j)+beta*B22(:,i,j),...
      DOFs,DOFs);
  end
end

%% JL20131016: TO BE REVISED FOR EFFECIENCY 
%% GlbMat: Term 2&3: Terms from integration by parts for Dirichlet edges 
A11 = zeros(NumDiriEgs,3,3);
B11 = zeros(NumDiriEgs,3,3);
GradBasFxn = [0,0; 1,0; 0,1];
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig=1:NumDiriEgs
  x1 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),1),1);
  y1 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),1),2);
  x2 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),2),1);
  y2 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),2),2);
  nml = repmat(TriMesh.NmlEg(DirichletEdge(ig),:),3,1);
  ie1 = TriMesh.edge2elem(DirichletEdge(ig),1);
  PermK1 = PermK(ie1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    X = qp(1);  Y = qp(2);
    BasFxn = [1; X; Y];
    Kdpgn1 = PermK1*dot(GradBasFxn,nml,2);
    for i=1:3
      for j=1:3
        A11(ig,i,j) = A11(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn(i)*Kdpgn1(j)*LenEg(DirichletEdge(ig));
      end
    end
  end
  % JL20131016: NEEDS DOUBLE-CHECK !!!
  B11(ig,:,:) = squeeze(A11(ig,:,:))';  
end
% Now assembling 
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
for i=1:3
  II1 = 3*(LblElt1-1) + i;
  for j=1:3
    JJ1 = 3*(LblElt1-1) + j;
    GlbMat = GlbMat + sparse(II1,JJ1,-A11(:,i,j)+beta*B11(:,i,j),...
      DOFs,DOFs);
  end
end

%% GlbMat: Term 4: Penalty on interior edges 
% JL20131016: TO BE REVISED FOR EFFECIENCY 
A11 = zeros(NumIntrEgs,3,3);  A12 = zeros(NumIntrEgs,3,3);  
A21 = zeros(NumIntrEgs,3,3);  A22 = zeros(NumIntrEgs,3,3);  
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig=1:NumIntrEgs
  x1 = TriMesh.node(TriMesh.edge(IntrEdge(ig),1),1);
  y1 = TriMesh.node(TriMesh.edge(IntrEdge(ig),1),2);
  x2 = TriMesh.node(TriMesh.edge(IntrEdge(ig),2),1);
  y2 = TriMesh.node(TriMesh.edge(IntrEdge(ig),2),2);
  % ie1 = TriMesh.edge2elem(IntrEdge(ig),1);
  % ie2 = TriMesh.edge2elem(IntrEdge(ig),2);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    X = qp(1);  Y = qp(2);
    BasFxn1 = [1; X; Y];
    BasFxn2 = [1; X; Y];
    for i=1:3
      for j=1:3
        A11(ig,i,j) = A11(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn1(i)*BasFxn1(j)*LenEg(IntrEdge(ig));
        A12(ig,i,j) = A12(ig,i,j)...
          - GAUSSQUAD.LINE(k,3)*BasFxn1(i)*BasFxn2(j)*LenEg(IntrEdge(ig));
        A21(ig,i,j) = A21(ig,i,j)...
          - GAUSSQUAD.LINE(k,3)*BasFxn2(i)*BasFxn1(j)*LenEg(IntrEdge(ig));
        A22(ig,i,j) = A22(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn2(i)*BasFxn2(j)*LenEg(IntrEdge(ig));
      end
    end
  end
end
LblElt1 = TriMesh.edge2elem(IntrEdge,1);
LblElt2 = TriMesh.edge2elem(IntrEdge,2);
for i=1:3
  II1 = 3*(LblElt1-1) + i;
  II2 = 3*(LblElt2-1) + i;
  for j=1:3
    JJ1 = 3*(LblElt1-1) + j;
    JJ2 = 3*(LblElt2-1) + j;
    GlbMat = GlbMat + sparse(II1,JJ1,A11(:,i,j)./LenEg(IntrEdge(ig))*alpha,...
      DOFs,DOFs);
    GlbMat = GlbMat + sparse(II1,JJ2,A12(:,i,j)./LenEg(IntrEdge(ig))*alpha,...
      DOFs,DOFs);
    GlbMat = GlbMat + sparse(II2,JJ1,A21(:,i,j)./LenEg(IntrEdge(ig))*alpha,...
      DOFs,DOFs);
    GlbMat = GlbMat + sparse(II2,JJ2,A22(:,i,j)./LenEg(IntrEdge(ig))*alpha,...
      DOFs,DOFs);
  end
end

%% JL20131016: TO BE REVISED FOR EFFECIENCY 
%% GlbMat: Term 5: Penalty on Dirichlet edges 
% For simplicity: a uniform penalty factor alpha 
A11 = zeros(NumDiriEgs,3,3);  
% Assuming NumDirichletEgs>0
for ig=1:NumDiriEgs
  x1 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),1),1);
  y1 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),1),2);
  x2 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),2),1);
  y2 = TriMesh.node(TriMesh.edge(DirichletEdge(ig),2),2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    X = qp(1);  Y = qp(2);
    BasFxn1 = [1; X; Y];
    for i=1:3
      for j=1:3
        A11(ig,i,j) = A11(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn1(i)*BasFxn1(j)*LenEg(DirichletEdge(ig));
      end
    end
  end
end
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
for i=1:3
  II1 = 3*(LblElt1-1) + i;
  for j=1:3
    JJ1 = 3*(LblElt1-1) + j;
    GlbMat = GlbMat + sparse(II1,JJ1,A11(:,i,j)./LenEg(DirichletEdge)*alpha,...
      DOFs,DOFs);
  end
end

%% Assembling the global right-hand side GlbRHS 

%% GlbRHS: Term 1: Integrals of source over triangular elements 
term1 = zeros(NumEms,3);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1)*TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2)*TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3)*TriMesh.node(TriMesh.elem(:,3),:);
  fval = EqnBC.fxnf(qp);
  X = qp(:,1);  Y = qp(:,2);
  term1(:,1) = term1(:,1) + GAUSSQUAD.TRIG(k,4)*fval.*area;
  term1(:,2) = term1(:,2) + GAUSSQUAD.TRIG(k,4)*(fval.*X).*area;
  term1(:,3) = term1(:,3) + GAUSSQUAD.TRIG(k,4)*(fval.*Y).*area;
end
GlbRHS = GlbRHS + term1;

%% GlbRHS: Term 2: Enforcing Neumman boundary condition 
if NumNeumEgs>0
  term2 = zeros(NumNeumEgs,3);
  x1 = TriMesh.node(TriMesh.edge(NeumannEdge,1),1);
  y1 = TriMesh.node(TriMesh.edge(NeumannEdge,1),2);
  x2 = TriMesh.node(TriMesh.edge(NeumannEdge,2),1);
  y2 = TriMesh.node(TriMesh.edge(NeumannEdge,2),2);
  LenNeumEg = sqrt((x2-x1).^2+(y2-y1).^2);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    uNval = EqnBC.fxnuN(qp);
    X = qp(:,1);  Y = qp(:,2);
    term2(:,1) = term2(:,1) + GAUSSQUAD.LINE(k,3)*uNval.*LenNeumEg;
    term2(:,2) = term2(:,2) + GAUSSQUAD.LINE(k,3)*uNval.*X.*LenNeumEg;
    term2(:,3) = term2(:,3) + GAUSSQUAD.LINE(k,3)*uNval.*Y.*LenNeumEg;
  end
  LblElt1 = TriMesh.edge2elem(NeumannEdge,1);
  for j=1:3 
    GlbRHS(LblElt1,j) = GlbRHS(LblElt1,j) - term2(:,j);  % Notice "-"
  end
end

%% GlbRHS: Term 3 & 4 could be merged ! 

%% GlbRHS: Term 3: Term from integration by parts for Dirichlet edges 
term3 = zeros(NumDiriEgs,3);
x1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),1);
x2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),1);
y1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),2);
y2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),2);
LenDiriEg = sqrt((x2-x1).^2+(y2-y1).^2);
nml = TriMesh.NmlEg(DirichletEdge,:);
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
PermK1 = PermK(LblElt1);
WO = ones(NumDiriEgs,1);  
ZO = zeros(NumDiriEgs,1);
GradBasFxn = zeros(NumDiriEgs,2,3);
Kdpgn = zeros(NumDiriEgs,3);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
  % X = qp(:,1);  Y = qp(:,2);
  GradBasFxn(:,:,2) = [WO,ZO];  
  GradBasFxn(:,:,3) = [ZO,WO];  
  pDval = EqnBC.fxnpD(qp);
  for j=1:3
    Kdpgn(:,j) = PermK1.*dot(squeeze(GradBasFxn(:,:,j)),nml,2);
    term3(:,j) = term3(:,j)...
      + GAUSSQUAD.LINE(k,3)*(Kdpgn(:,j).*pDval).*LenDiriEg;
  end
end
for j=1:3 
  GlbRHS(LblElt1,j) = GlbRHS(LblElt1,j) + term3(:,j)*beta;  % merge?
end

%% GlbRHS: Term 4: Penalty on Dirichlet edges 
% For simplicity: a uniform penalty factor alpha 
% JL20131016: TO BE REVISED FOR EFFECIENCY 
% Assuming NumDirichletEgs>0
term4 = zeros(NumDiriEgs,3);
x1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),1);
y1 = TriMesh.node(TriMesh.edge(DirichletEdge,1),2);
x2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),1);
y2 = TriMesh.node(TriMesh.edge(DirichletEdge,2),2);
LenDiriEg = sqrt((x2-x1).^2+(y2-y1).^2);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
  pDval = EqnBC.fxnpD(qp);
  X = qp(:,1);  Y = qp(:,2);
  term4(:,1) = term4(:,1) + GAUSSQUAD.LINE(k,3)*pDval.*LenDiriEg;
  term4(:,2) = term4(:,2) + GAUSSQUAD.LINE(k,3)*(pDval.*X).*LenDiriEg;
  term4(:,3) = term4(:,3) + GAUSSQUAD.LINE(k,3)*(pDval.*Y).*LenDiriEg;
end
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
for j=1:3 
  GlbRHS(LblElt1,j) = GlbRHS(LblElt1,j) + term4(:,j)./LenDiriEg*alpha;
end

%% Solving the global lin. sys. directly 
GlbRHS = reshape(GlbRHS',3,NumEms);
GlbRHS = GlbRHS(:);
sln = GlbMat\GlbRHS;
% sln = zeros(DOFs,1);

%% Converting 
CofNS = (reshape(sln,3,NumEms))';

return;