function [Q,G,H,R]=assemb(varargin) % bl,p,e,u,time,sdl
%ASSEMB Assembles boundary condition contributions in a PDE problem.
%
%       [Q,G,H,R]=ASSEMB(B,P,E) assembles the matrices Q and
%       H, and the vectors G and R.
%       Q should be added to the system matrix and contains contributions
%       from Mixed boundary conditions.
%       G should be added to the right-hand side and contains contributions
%       from Neumann and Mixed boundary conditions.
%       H*u=R represents the Dirichlet type boundary conditions.
%
%       The following call sequences are also allowed:
%       [Q,G,H,R]=ASSEMB(B,P,E,U0)
%       [Q,G,H,R]=ASSEMB(B,P,E,U0,TIME)
%       [Q,G,H,R]=ASSEMB(B,P,E,U0,TIME,SDL)
%       [Q,G,H,R]=ASSEMB(B,P,E,TIME)
%       [Q,G,H,R]=ASSEMB(B,P,E,TIME,SDL)
%
%       [Q,G,H,R]=ASSEMB(PDEM,...) This format allows the boundary conditions 
%       B and the mesh nodes and associativity, P and E to be defined by 
%       PDEM a PDEMODEL object that contains B, P and E.
%
%       The input parameters B, P, E, U0, TIME, and SDL have the same
%       meaning as in ASSEMPDE.
%
%       See also CREATEPDE, pde.PDEMODEL, ASSEMPDE, PDEBOUND

%       Copyright 1994-2020 The MathWorks, Inc.

if nargin > 0
    [varargin{:}] = convertStringsToChars(varargin{:});
end

narginchk(1,6);
gotu=false;
gottime=false;
gotsdl=false;

narg = nargin;
offset=0;
meshedpde=false;

if isa(varargin{1}, 'pde.PDEModel') 
    mypde = varargin{1}; 
    if isa(mypde.Mesh, 'pde.FEMesh')
      meshedpde=true;
    elseif ~mypde.IsTwoD
      error(message('pde:assemb:pdemodelNoMesh'));  
    end
    if nargin==2 && ~isfloat(varargin{2})
       error(message('pde:assemb:nonNumericTimeU0')); 
    end
end
if meshedpde
    mypde = varargin{1};
    msh = mypde.Mesh;
    [p,e,~] = msh.meshToPet();
    bl = pde.PDEModel(mypde.PDESystemSize);
    if isempty(mypde.BoundaryConditions)
        bc = [];
    else
        bc = mypde.BoundaryConditions.BoundaryConditionAssignments;
    end
    offset=2;
    narg = narg+2;
else
    narginchk(3,6);
    bl = varargin{1};
    if isa(bl, 'pde.PDEModel')
        if isempty(bl.BoundaryConditions)
            bc = [];
        else
            bc = bl.BoundaryConditions.BoundaryConditionAssignments;
        end
    end
    p = varargin{2};
    e = varargin{3};
end


if narg==3
  % No action
elseif narg==4
  if size(varargin{4-offset},1)>1
    u = varargin{4-offset};
    gotu=true;
  else
    time=varargin{4-offset};
    gottime=true;
  end
elseif narg==5
  if size(varargin{4-offset},1)>1
    gotu=true;
    gottime=true;
    u = varargin{4-offset};
    time = varargin{5-offset};
  else
    sdl=varargin{5-offset};
    time=varargin{4-offset};
    gottime=true;
    gotsdl=true;
  end
elseif narg==6
  u = varargin{4-offset};
  time = varargin{5-offset};
  sdl = varargin{6-offset};
  gotu=true;
  gottime=true;
  gotsdl=true;
else
  error(message('pde:assema:nargin'));
end

if( isa(bl, 'pde.PDEModel'))
    npr = size(p,1);
    if ~gotu
        u = [];
    end
    if ~gottime
        time = [];
    end
    if(npr == 3)
        if gotsdl && ~isempty(sdl)
            error(message('pde:pdeModel:sdl3D'));
        end
    end
    bcImpl = pde.internal.pde3DBCImpl(bl,bc,false,true,p,e); % false for OLD BC implementation
    bcmat = bcImpl.getBCMatrices(u,time,'HRQG');
    Q = bcmat.Q;
    G = bcmat.G;
    H = bcmat.H;
    R = bcmat.R;
    return;
end

pdeIsValid.e(e); % check the validity of the e-matrix
if ~gotsdl
  ie=pdesde(e);
else
  ie=pdesde(e,sdl);
end

np=size(p,2); % Number of points

if gotu
  N=size(u,1)/np;
end

% Try to calculate N by evaluating the boundary conditions on the first boundary edge
if isempty(ie) && ~gotu
  % Look at the first edge
  [q,g,h,r] = pdeexpd(p,e(:,1),bl);
  N = size(g,1);
end

% Nothing to do?
if isempty(ie)
  Q=sparse(N*np,N*np);
  G=sparse(N*np,1);
  H=sparse(N*np,N*np);
  R=sparse(N*np,1);
  % Done
  return;
end

% Get rid of unwanted edges
e=e(:,ie);

ne=size(e,2); % Number of edges

if ~gotu && ~gottime
  [q,g,h,r]=pdeexpd(p,e,bl);
elseif gotu && ~gottime
  [q,g,h,r]=pdeexpd(p,e,u,bl);
elseif ~gotu && gottime
  [q,g,h,r]=pdeexpd(p,e,time,bl);
else % gotu & gottime
  [q,g,h,r]=pdeexpd(p,e,u,time,bl);
end

if ~gotu
  N=size(g,1);
end

Q=sparse(N*np,N*np);
G=sparse(N*np,1);
H=sparse(N*np,N*np);
R=sparse(N*np,1);

% ie1 and ie2 contain the beginning and ending node numbers for each
% edge on the boundary, respectively
ie1=e(1,:);
ie2=e(2,:);
d=sqrt((p(1,ie1)-p(1,ie2)).^2+(p(2,ie1)-p(2,ie2)).^2); % Side length

% G
if any(size(g)-[N ne])
  error(message('pde:assemb:SizeG'));
end
if any(g(:)) | any(isnan(g(:)))
  g=(g.*(ones(N,1)*d/2));
  for k=1:N
	G=G+sparse(ie1+(k-1)*np,1,g(k,:),N*np,1);
    G=G+sparse(ie2+(k-1)*np,1,g(k,:),N*np,1);
  end
end

% Q
if any(size(q)-[N^2 ne])
  error(message('pde:assemb:SizeQ'));
end
if any(q(:)) | any(isnan(q(:)))
  qd=(q.*(ones(N*N,1)*d/3)); % Diagonal values
  qod=0.5*qd; % Off diagonal values
  m=1;
  for l=1:N
    for k=1:N
      Q=Q+sparse(ie1+(k-1)*np,ie2+(l-1)*np,qod(m,:),N*np,N*np);
      Q=Q+sparse(ie2+(k-1)*np,ie1+(l-1)*np,qod(m,:),N*np,N*np);
      Q=Q+sparse(ie1+(k-1)*np,ie1+(l-1)*np,qd(m,:),N*np,N*np);
      Q=Q+sparse(ie2+(k-1)*np,ie2+(l-1)*np,qd(m,:),N*np,N*np);
	  m=m+1;
    end
  end
end

if any(size(h)-[N^2 2*ne])
  error(message('pde:assemb:SizeH'));
end
if any(size(r)-[N 2*ne])
  error(message('pde:assemb:SizeR'));
end
if any(h(:)) || any(isnan(h(:)))
  % the h matrix has two columns for each boundary edge in the FE
  % model. The first ne columns contain the input h-coefficients at
  % the first node of each edge and the second ne columns contain the
  % h-coefficients at the second node of each edge.
  
  % Preprocess h & r

  % The first ne entries in ie are the node numbers of the first node
  % in each FE boundary edge. The next ne entries are node numbers of
  % the second node in each edge.
  ie=[ie1 ie2];
  % Order h and r by node number
  [ie,k]=sort(ie);
  h=h(:,k);
  r=r(:,k);

  % Now remove redundant entries
  
  % a BC is duplicate only if all the h-coefficients are the same and
  % the node numbers are the same
  % Have to handle the case with NaN specially. This is used to see if
  % an expression is function of some variable (e.g. time).
  isNanH = isnan(h);
  hCompare = h;
  hCompare(isNanH) = 1.0;
  dBcDef = [hCompare; ie];
  xtmp = diff(dBcDef, 1, 2);
  k = sum(abs(xtmp));
  kk=find(k);
  k(kk)=ones(size(kk));
  iie=cumsum([1 k]);
  aa=sparse(1:length(ie),iie,0.5);
  ie=ie([kk length(ie)]);
  h=h*aa;
  r=r*aa;

  nie=length(ie);

  % Spread out to one equation per column
  ie=ones(N,1)*ie;
  ie=ie(:)';
  i1=(1:N)'*ones(1,N);
  i1=ones(N*N,1)*N*(0:nie-1)+i1(:)*ones(1,nie);
  i2=ones(N,1)*(1:N);
  i2=i2(:)*ones(1,nie);
  hh=zeros(N,N*nie);
  hh(i2+N*(i1-1))=h(:);
  h=hh;
  r=r(:).';
  nie=length(ie);

  % Remove more redundancy
  if N==1
    k=find(h);
  else
    k=find(max(abs(h)));
  end
  ie=ie(k);
  h=h(:,k);
  r=r(:,k);
  nie=length(ie);

  % H
  i1=ones(N,1)*(1:nie);
  i2=ones(N,1)*ie+(0:N-1)'*np*ones(1,nie);
  H=sparse(i1,i2,h,nie,N*np);

  % R
  
  R=sparse(1:nie,1,r,nie,1);

end

