function [sln] = LinElas_CG_TriEP1_AsmSlv( ...
  EqnBC,BndryCondType,TriMesh,GAUSSQUAD)
%% Solving linear elasticity by CG EP1 elements on a triangular mesh
% EP1: CP1^2 plus 3 edge bubble functions (using edge normal vectors) 
% See 
% See also Bernardi,Ragel, MCOM, 44(1985), pp.71-79. 
% James Liu, ColoState; 2012/07--2020/10 

%% Equation info 
lambda = EqnBC.lambda;
    mu = EqnBC.mu;

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

%% Sorting out boundary edges: 1=Dirichlet, 2=Neumann 
DirichletEdge = find(BndryCondType(TriMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(TriMesh.BndryEdge+1)==2);

%% Boundary info 
NumDirichletEgs = size(DirichletEdge,1);
NumNeumannEgs = size(NeumannEdge,1);
% JL20160808: TO BE REVISED FOR EFFICIENCY !
DirichletNodeFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgs
  j1 = TriMesh.edge(DirichletEdge(ig),1);
  j2 = TriMesh.edge(DirichletEdge(ig),2);
  DirichletNodeFlag(j1) = 1;
  DirichletNodeFlag(j2) = 1;
end
DirichletNode = find(DirichletNodeFlag);

%% Setup for discrete linear system 
DOFs = 2*NumNds + NumEgs;
% GlbMat = sparse(DOFs,DOFs);
% GlbRHS = zeros(DOFs,1);

%% Computing element matrices 
% Element matrix: strain-strain 
EltMatStrnStrn = zeros(NumEms,9,9);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  StrnPhi = zeros(NumEms,9,2,2);
  StrnPhi(:,:,1,1) = BFGrad(:,:,1,1);
  StrnPhi(:,:,2,2) = BFGrad(:,:,2,2);
  StrnPhi(:,:,1,2) = 0.5*(BFGrad(:,:,1,2)+BFGrad(:,:,2,1));
  StrnPhi(:,:,2,1) = StrnPhi(:,:,1,2);
  for i=1:9
    for j=1:9
      EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j)...
        + ( StrnPhi(:,i,1,1).*StrnPhi(:,j,1,1)...
          + StrnPhi(:,i,1,2).*StrnPhi(:,j,1,2)...
          + StrnPhi(:,i,2,1).*StrnPhi(:,j,2,1)...
          + StrnPhi(:,i,2,2).*StrnPhi(:,j,2,2)) .* area * GAUSSQUAD.TRIG(k,4);
    end
  end
end
% Element matrix: div(avg)-div(avg) 
EltMatDivDiv = zeros(NumEms,9,9);
EltDivAvg = zeros(NumEms,9);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
for j=1:3
  EltDivAvg(:,2*j-1) = CP1BasGrad(:,j,1);
  EltDivAvg(:,2*j  ) = CP1BasGrad(:,j,2);
end
EltDivAvg(:,7) = TriMesh.LenEg(TriMesh.elem2edge(:,1))./(6*area);
EltDivAvg(:,8) = TriMesh.LenEg(TriMesh.elem2edge(:,2))./(6*area);
EltDivAvg(:,9) = TriMesh.LenEg(TriMesh.elem2edge(:,3))./(6*area);
for i=1:9
  for j=1:9
    EltMatDivDiv(:,i,j) = EltDivAvg(:,i) .* EltDivAvg(:,j) .* area;
  end
end
% Altogether 
% EltMat = zeros(NumEms,9,9);
EltMat = (2*mu)*EltMatStrnStrn + lambda*EltMatDivDiv;
% Modification with signs for edge normal vectors 
for i=1:3
  for j=1:9
    EltMat(:,6+i,j) = EltMat(:,6+i,j) .* TriMesh.EmEgSign(:,i);
  end
end
for i=1:9
  for j=1:3
    EltMat(:,i,6+j) = EltMat(:,i,6+j) .* TriMesh.EmEgSign(:,j);
  end
end

%% Assembling GlbMat 
GlbMat = sparse(DOFs,DOFs);
pos = zeros(NumEms,9);
for j=1:3  % Treating vertices and edges in one loop 
  k = TriMesh.elem(:,j);
  pos(:,2*j-1) = 2*k - 1;
  pos(:,2*j  ) = 2*k;
  k = TriMesh.elem2edge(:,j);
  pos(:,6+j) = 2*NumNds + k;
end
% Now assembling 
for i=1:9 
  II = pos(:,i);
  for j=1:9 
    JJ = pos(:,j);
    GlbMat = GlbMat + sparse(II,JJ, EltMat(:,i,j), DOFs,DOFs);
  end
end

%% Assembling GlbRHS 
GlbRHS = zeros(DOFs,1);
% Computing elementwise contributions of 9 basis functions 
EltRHS = zeros(NumEms,9);
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),:);
  BFVal = CG_TriBR1_BasFxnVal(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  fval = EqnBC.fxnf(qp);
  for j=1:9
    EltRHS(:,j) = EltRHS(:,j) + GAUSSQUAD.TRIG(k,4)...
      *( BFVal(:,j,1).*fval(:,1) + BFVal(:,j,2).*fval(:,2) ).*area;
  end  
end
% Modification for 3 edge-based basis functions 
for j=1:3
  EltRHS(:,6+j) = EltRHS(:,6+j) .* TriMesh.EmEgSign(:,j);
end
% Assembling contributions of elementwise 9 basis functions 
for ie=1:NumEms
  for j=1:3  % Treating vertices and edges in one loop 
    k = TriMesh.elem(ie,j);
    GlbRHS(2*k-1) = GlbRHS(2*k-1) + EltRHS(ie,2*j-1);
    GlbRHS(2*k  ) = GlbRHS(2*k  ) + EltRHS(ie,2*j  );
    k = TriMesh.elem2edge(ie,j);
    GlbRHS(2*NumNds+k) = GlbRHS(2*NumNds+k) + EltRHS(ie,6+j);
  end
end

%% JL20171221: TO BE FINISHED/REVISED BY Graham 
%% Incorporating boundary conditions: Neumann as natural 
if NumNeumannEgs>0
  % For end-nodes 
  NeumannBC = zeros(NumNeumannEgs,2,2);
  j1 = TriMesh.edge(NeumannEdge,1);
  j2 = TriMesh.edge(NeumannEdge,2);
  x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2);
  x2 = TriMesh.node(j2,1);  y2 = TriMesh.node(j2,2);
  LenEgNeumann = 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];
    tNval = EqnBC.fxntN(qp);
    NeumannBC(:,1,1) = NeumannBC(:,1,1) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,1) * tNval(:,1) .* LenEgNeumann;
    NeumannBC(:,1,2) = NeumannBC(:,1,2) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,1) * tNval(:,2) .* LenEgNeumann;
    NeumannBC(:,2,1) = NeumannBC(:,2,1) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,2) * tNval(:,1) .* LenEgNeumann;
    NeumannBC(:,2,2) = NeumannBC(:,2,2) + GAUSSQUAD.LINE(k,3)...
      * GAUSSQUAD.LINE(k,2) * tNval(:,2) .* LenEgNeumann;
  end
  for ig=1:NumNeumannEgs 
    j1 = TriMesh.edge(NeumannEdge(ig),1);
    j2 = TriMesh.edge(NeumannEdge(ig),2);
    GlbRHS(2*j1-1) = GlbRHS(2*j1-1) - NeumannBC(ig,1,1);
    GlbRHS(2*j1  ) = GlbRHS(2*j1  ) - NeumannBC(ig,1,2);
    GlbRHS(2*j2-1) = GlbRHS(2*j2-1) - NeumannBC(ig,2,1);
    GlbRHS(2*j2  ) = GlbRHS(2*j2  ) - NeumannBC(ig,2,2);
  end
  % For edges 
  NeumannBC = zeros(NumNeumannEgs,1);
  j1 = TriMesh.edge(NeumannEdge,1);
  j2 = TriMesh.edge(NeumannEdge,2);
  x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2);
  x2 = TriMesh.node(j2,1);  y2 = TriMesh.node(j2,2);
  LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
  nml = TriMesh.EgNml(NeumannEdge,:);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    tNval = EqnBC.fxntN(qp);
    psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
    NeumannBC = NeumannBC + GAUSSQUAD.LINE(k,3)...
      * sum(tNval.*nml,2) * psi .* LenEgNeumann;
  end
  GlbRHS(2*NumNds+NeumannEdge) = GlbRHS(2*NumNds+NeumannEdge) - NeumannBC;
end

%% Incorporating boundary conditions: Dirichlet as essential 
sln = zeros(DOFs,1);
% Assuming NumDirichletEgs>0 
% For nodes 
uDval = EqnBC.fxnuD(TriMesh.node(DirichletNode,:));
sln(2*DirichletNode-1) = uDval(:,1);
sln(2*DirichletNode  ) = uDval(:,2);
% For edges 
intgrl1 = zeros(NumDirichletEgs,1);
intgrl2 = zeros(NumDirichletEgs,1);
% intplt = zeros(NumDirichletEgs,1);
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);
uD1 = EqnBC.fxnuD([x1,y1]); 
uD2 = EqnBC.fxnuD([x2,y2]); 
% LenDirichletEdge = sqrt((x2-x1).^2+(y2-y1).^2);
nml = TriMesh.EgNml(DirichletEdge,:);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts 
  qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
  uD = EqnBC.fxnu(qp);
  psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
  intplt = GAUSSQUAD.LINE(k,1)*uD1 + GAUSSQUAD.LINE(k,2)*uD2;
  res = uD - intplt;
  intgrl1 = intgrl1 + GAUSSQUAD.LINE(k,3) * psi;
  intgrl2 = intgrl2 + GAUSSQUAD.LINE(k,3) * sum(res.*nml,2);
end
sln(2*NumNds+DirichletEdge) = intgrl2./intgrl1;

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

%% Flag 
flag = ones(DOFs,1);
flag(2*DirichletNode-1) = 0;
flag(2*DirichletNode  ) = 0;
flag(2*NumNds+DirichletEdge) = 0;
FreeNdEg = find(flag);

%% Solving the reduced linear system directly 
sln(FreeNdEg) = GlbMat(FreeNdEg,FreeNdEg) \ GlbRHS(FreeNdEg);

return;