function [G,H,R]=lassemb(b,p,e)
% function [G,H,R]=lassemb(b,p,e);
%
% Handles the contributions of boundary conditions given in b.
% Assumes h=1 and q=0, that is, conditions of type
%   u = r      Dirichlet  or
%   u' = g     Neumann.
% No reduction to the size of the system is done.

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

G=sparse(np,1);
R=sparse(np,1);
H=sparse(np,np);

ie=pdesde(e);  % Indices of outer boundary edges see help pdesde
% Nothing to do?
if length(ie)==0
  % Done
  return;
end

% Get rid of unwanted edges
e=e(:,ie);
ne=size(e,2);  % Number of edges

%%% Diriclet conditions
%%%%%%%%%%%%%%%%%%%%%%%

id=find(b(2,:)==1);  % Dirichlet cond. on these boundary segments
for i=id
  ec = find(e(5,:)==i);     % edges of current boundary segment
  pc = e(1:2,ec); pc=pc(:); % points of current boundary segment
  npc= size(pc,1);          % number of points
  
  x=p(1,pc); 
  y=p(2,pc);
  str=char(b(10:10+b(6,i)-1,i));  % See this, helps in using b-matrix
  r=eval(str);                    % u = r
  if ~(size(r,2) == npc)          % if b.c. is constant: r='1'
    r=r*ones(npc,1);
  end
  for j=1:npc
    R(pc(j))=R(pc(j))+r(j)/2;
    H(pc(j),pc(j))=H(pc(j),pc(j))+0.5;
  end
  
end;

% reshaping the matrices H and R
[i j]=find(H);                 
ndp=size(i,1);                 % number of Dirichlet points
h=sparse(ndp,np);
h(1:ndp,j)=H(i,j); H=h;
r=sparse(ndp,1);
r=R(i); R=r;                   

%%% Neumann conditions
%%%%%%%%%%%%%%%%%%%%%%

in=find(b(2,:)==0); % Neumann cond. on these boundary segments
for i=in
  ec = find(e(5,:)==i);  % edges of current boundary segment
  for j=ec
    pc = e(1:2,j);       % points of current boundary segment
    npc= 2;              % number of points

    xp=p(1,pc); yp=p(2,pc);
    x=mean(xp); y=mean(yp);        % mid-point Pb of the edge    

    lg=b(4,i);                     % length of string reprsenting g
    lq=b(3,i);
    str=char(b(5+lq:5+lq+lg-1,i)); % u' = g  
    g=eval(str);                   % g(Pb)
    
    Gt = g/2*norm([xp(1) yp(1)]-[xp(2) yp(2)])*ones(2,1);
    G(pc)=G(pc)+Gt;
  end
end


%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-05-13 15:39:37 EEST

