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

%% Mesh info 
NumEms = TriMesh.NumEms;
% NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag>=2
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 = NumEms*6;
GlbMat = sparse(DOFs,DOFs);
GlbRHS = zeros(NumEms,6);  % This is nonconventional! 

%% Assembling the global coeff. matrix (GlbMat) 

%% GlbMat: Term 1: grad-grad on all triangular elements 
ZO = zeros(NumEms,1);
IX = ZO;  IY = ZO;  IX2 = ZO;  IXY = ZO;  IY2 = ZO;
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),:);
  X = qp(:,1);  Y = qp(:,2);
  IX = IX + GAUSSQUAD.TRIG(k,4)*X.*area;
  IY = IY + GAUSSQUAD.TRIG(k,4)*Y.*area;
  IX2 = IX2 + GAUSSQUAD.TRIG(k,4)*(X.^2).*area;
  IXY = IXY + GAUSSQUAD.TRIG(k,4)*(X.*Y).*area;
  IY2 = IY2 + GAUSSQUAD.TRIG(k,4)*(Y.^2).*area;
end
MatGG = zeros(NumEms,6,6);
for j=1:6; MatGG(:,1,j)=ZO; end
MatGG(:,2,2) = area;  MatGG(:,2,4) = 2*IX;  MatGG(:,2,5) = IY;
MatGG(:,3,3) = area;  MatGG(:,3,5) = IX;  MatGG(:,3,6) = 2*IY;
MatGG(:,4,4) = 4*IX2;  MatGG(:,4,5) = 2*IXY;
MatGG(:,5,5) = IX2+IY2;  MatGG(:,5,6) = 2*IXY;  
MatGG(:,6,6) = 4*IY2;
for i=2:6
  for j=1:(i-1)
    MatGG(:,i,j) = MatGG(:,j,i);
  end
end
% Now assembling 
% Assuming PermK is an elementwise scalar constant 
for i=1:6
  II = 6*(0:NumEms-1) + i;
  for j=1:6
    JJ = 6*(0:NumEms-1) + j;
    Aij = PermK.*MatGG(:,i,j);
    GlbMat = GlbMat + sparse(II,JJ,Aij,DOFs,DOFs);
  end
end
clear MatGG Aij;

%% GlbMat: Term 2&3: Terms from integration by parts for interior edges 
% JL20131010: TO BE REVISED FOR EFFECIENCY 
A11 = zeros(NumIntrEgs,6,6);  A12 = zeros(NumIntrEgs,6,6);  
A21 = zeros(NumIntrEgs,6,6);  A22 = zeros(NumIntrEgs,6,6);  
B11 = zeros(NumIntrEgs,6,6);  B12 = zeros(NumIntrEgs,6,6);  
B21 = zeros(NumIntrEgs,6,6);  B22 = zeros(NumIntrEgs,6,6);  
% BasFxn = zeros(6,1);
% GradBasFxn = zeros(6,2);
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),:),6,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; X^2; X*Y; Y^2];
    GradBasFxn = [0,0; 1,0; 0,1; 2*X,0; Y,X; 0,2*Y];
    Kdpgn1 = PermK1*dot(GradBasFxn,nml,2);
    Kdpgn2 = PermK2*dot(GradBasFxn,nml,2);
    for i=1:6
      for j=1:6
        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
  % JL20131010: 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:6
  II1 = 6*(LblElt1-1) + i;
  II2 = 6*(LblElt2-1) + i;
  for j=1:6
    JJ1 = 6*(LblElt1-1) + j;
    JJ2 = 6*(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

%% GlbMat: Term 2&3: Terms from integration by parts for Dirichlet edges 
% JL20131010: TO BE REVISED FOR EFFECIENCY 
A11 = zeros(NumDiriEgs,6,6);
B11 = zeros(NumDiriEgs,6,6);
% BasFxn = zeros(6,1);
% GradBasFxn = zeros(6,2);
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),:),6,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; X^2; X*Y; Y^2];
    GradBasFxn = [0,0; 1,0; 0,1; 2*X,0; Y,X; 0,2*Y];
    Kdpgn1 = PermK1*dot(GradBasFxn,nml,2);
    for i=1:6
      for j=1:6
        A11(ig,i,j) = A11(ig,i,j)...
          + GAUSSQUAD.LINE(k,3)*BasFxn(i)*Kdpgn1(j)*LenEg(DirichletEdge(ig));
      end
    end
  end
  % JL20131010: NEEDS DOUBLE-CHECK !!!
  B11(ig,:,:) = squeeze(A11(ig,:,:))';  
end
% Now assembling 
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
for i=1:6
  II1 = 6*(LblElt1-1) + i;
  for j=1:6
    JJ1 = 6*(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 
% JL20131017: TO BE REVISED FOR EFFECIENCY 
A11 = zeros(NumIntrEgs,6,6);  A12 = zeros(NumIntrEgs,6,6);  
A21 = zeros(NumIntrEgs,6,6);  A22 = zeros(NumIntrEgs,6,6);  
% BasFxn1 = zeros(6,1);  % 1,x,y,x^2,xy,y^2 
% BasFxn2 = zeros(6,1);  % 1,x,y,x^2,xy,y^2 
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);
  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; X^2; X*Y; Y^2];
    BasFxn2 = [1; X; Y; X^2; X*Y; Y^2];
    for i=1:6
      for j=1:6
        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:6
  II1 = 6*(LblElt1-1) + i;
  II2 = 6*(LblElt2-1) + i;
  for j=1:6
    JJ1 = 6*(LblElt1-1) + j;
    JJ2 = 6*(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

%% GlbMat: Term 5: Penalty on Dirichlet edges  
% For simplicity: a uniform penalty factor alpha 
% JL20131016: TO BE REVISED FOR EFFECIENCY 
A11 = zeros(NumDiriEgs,6,6);  
% 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; X^2; X*Y; Y^2];
    for i=1:6
      for j=1:6
        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:6
  II1 = 6*(LblElt1-1) + i;
  for j=1:6
    JJ1 = 6*(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,6);
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;
  term1(:,4) = term1(:,4) + GAUSSQUAD.TRIG(k,4)*(fval.*(X.^2)).*area;
  term1(:,5) = term1(:,5) + GAUSSQUAD.TRIG(k,4)*(fval.*(X.*Y)).*area;
  term1(:,6) = term1(:,6) + GAUSSQUAD.TRIG(k,4)*(fval.*(Y.^2)).*area;
end
GlbRHS = GlbRHS + term1;

%% GlbRHS: Term 2: Enforcing Neumman boundary condition 
if NumNeumEgs>0
  term2 = zeros(NumNeumEgs,6);
  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;
    term2(:,4) = term2(:,4) + GAUSSQUAD.LINE(k,3)*uNval.*(X.^2).*LenNeumEg;
    term2(:,5) = term2(:,5) + GAUSSQUAD.LINE(k,3)*uNval.*(X.*Y).*LenNeumEg;
    term2(:,6) = term2(:,6) + GAUSSQUAD.LINE(k,3)*uNval.*(Y.^2).*LenNeumEg;
  end
  LblElt1 = TriMesh.edge2elem(NeumannEdge,1);
  for j=1:6 
    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,6);
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);
nml = TriMesh.NmlEg(DirichletEdge,:);
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
PermK1 = PermK(LblElt1);
WO = ones(NumDiriEgs,1);  
ZO = zeros(NumDiriEgs,1);
GradBasFxn = zeros(NumDiriEgs,2,6);
Kdpgn = zeros(NumDiriEgs,6);
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(:,:,1) = [ZO,ZO];  GradBasFxn(:,:,2) = [WO,ZO]; 
  GradBasFxn(:,:,3) = [ZO,WO];  GradBasFxn(:,:,4) = [2*X,ZO]; 
  GradBasFxn(:,:,5) = [Y,X];    GradBasFxn(:,:,6) = [ZO,2*Y];
  pDval = EqnBC.fxnpD(qp);
  for j=1:6
    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:6 
  GlbRHS(LblElt1,j) = GlbRHS(LblElt1,j) + beta*term3(:,j);
end

%% GlbRHS: Term 4: Penalty on Dirichlet edges 
% For simplicity: a uniform penalty factor alpha 
% JL20131009: TO BE REVISED FOR EFFECIENCY 
% Assuming NumDirichletEgs>0
term4 = zeros(NumDiriEgs,6);
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;
  term4(:,4) = term4(:,4) + GAUSSQUAD.LINE(k,3)*(pDval.*(X.^2)).*LenDiriEg;
  term4(:,5) = term4(:,5) + GAUSSQUAD.LINE(k,3)*(pDval.*(X.*Y)).*LenDiriEg;
  term4(:,6) = term4(:,6) + GAUSSQUAD.LINE(k,3)*(pDval.*(Y.^2)).*LenDiriEg;
end
LblElt1 = TriMesh.edge2elem(DirichletEdge,1);
for j=1:6 
  GlbRHS(LblElt1,j) = GlbRHS(LblElt1,j) + term4(:,j)./LenDiriEg*alpha;
end

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

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

return;