function [Q,G,H,R]=assemb(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)
%
%       The input parameters B, P, E, U0, TIME, and SDL have the same
%       meaning as in ASSEMPDE.
%
%       See also ASSEMPDE, PDEBOUND

%       A. Nordmark 11-08-94, AN 12-20-94.
%       Copyright (c) 1994-98 by The MathWorks, Inc.
%       $Revision: 1.7 $  $Date: 1997/11/25 17:04:40 $


gotu=0;
gottime=0;
gotsdl=0;

if nargin==3
  % No action
elseif nargin==4
  if size(u,1)>1
    gotu=1;
  else
    time=u;
    gottime=1;
  end
elseif nargin==5
  if size(u,1)>1
    gotu=1;
    gottime=1;
  else
    sdl=time;
    time=u;
    gottime=1;
    gotsdl=1;
  end
elseif nargin==6,
  gotu=1;
  gottime=1;
  gotsdl=1;
else
  error('Wrong number of input arguments');
end

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

% Stupid case where we must determine the number of variables in a kludgy way
if length(ie)==0 & ~gotu
  % Look at the first edge
  [q,g,h,r]=pdeexpd(p,e(:,1),bl);
  N=size(g,1);
end

% Nothing to do?
if length(ie)==0
  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

avoid_arithmetics_on_NaNs=0;
if gottime
  if ~isieee & any(isnan(time))
    % PARABOLIC and HYPERBOLIC probes with time=NaN for
    % explicit time dependence. We must not cause any floating point
    % exceptions in this case.
   avoid_arithmetics_on_NaNs=1;
  end
end

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

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('g has wrong size');
end
if any(g(:))
  if avoid_arithmetics_on_NaNs & any(isnan(g))
    G=sparse(1,1,NaN,N*np,1);
  else
    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
end

% Q
if any(size(q)-[N^2 ne])
  error('q has wrong size');
end
if any(q(:))
  if avoid_arithmetics_on_NaNs & any(isnan(q))
    Q=sparse(1,1,NaN,N*np,N*np);
  else
    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
end

if any(size(h)-[N^2 2*ne])
  error('h has wrong size');
end
if any(size(r)-[N 2*ne])
  error('r has wrong size');
end
if any(h(:))

  % Preprocess h & r

  ie=[ie1 ie2];
  % Now remove redundant entries
  [ie,k]=sort(ie);
  h=h(:,k);
  r=r(:,k);

  k=diff(ie);
  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)]);
  if avoid_arithmetics_on_NaNs & any(isnan(h))
    h=NaN*(h~=0);
    h=h(:,[kk length(ie)]);
  else
    h=h*aa;
  end
  if avoid_arithmetics_on_NaNs & any(isnan(r))
    r=NaN*(r~=0);
    r=r(:,[kk length(ie)]);
  else
    r=r*aa;
  end

  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
  if avoid_arithmetics_on_NaNs & any(isnan(h))
    H=sparse(1,1,NaN,nie,N*np);
  else
    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);
  end

  % R
  if avoid_arithmetics_on_NaNs & any(isnan(r))
    R=sparse(1,1,NaN,nie,1);
  else
    R=sparse(1:nie,1,r,nie,1);
  end

end


