% BEMLAP-MAT project
% Matlab/Freemat codes
% Copyright 2008 Stephen Kirkup 
% Issued under the GNU General Public License 2007, see gpl.txt
% www.kirkup.info/opensource
% www.boundary-element-method.com

% function [phi_D,phi_S,v_S] = libem2_indirect(n_D,p_D,n_S,vertpts_S,elemvert_S,alpha_S,beta_S,f_S)
% Returns the solution phi_D at the n_D domain points p_D and at the collocation points.  
% The boundary is made up of n elements; vertpts lists the coordinates
%  of the edges of the elements and elemvert lists the pairs of vertices
%  that define each element. alpha_S, beta_S and f_S determine the Robin boundary condition.

function [phi_D,phi_S,v_S] =libem2_indirect(n_D,p_D,n_S,vertpts_S,elemvert_S,alpha_S,beta_S,f_S)
   
% calculate the boundary density function sigma

%  calculate L_SS, M_SS, Mt_SS and N_SS, the discrete form of the operators for the collocation points
   [L_SS,M_SS,Mt_SS,N_SS] =lbem2_on(n_S,vertpts_S,elemvert_S,true,true,true,true);
   Mt_SSplus=Mt_SS+eye(n_S)/2;
   M_SSminus=M_SS-eye(n_S)/2;
   mu=norm(Mt_SSplus)/norm(N_SS);
   matrix1=L_SS+mu*M_SSminus;
   matrix2=Mt_SSplus+mu*N_SS;
   for (i=1:n_S)
      matrix(i,:)=alpha_S(i)*matrix1(i,:)+beta_S(i)*matrix2(i,:);
   end
   sigma=matrix\f_S';
 
 
% calculate the solution on the boundary (often not necessary)
   phi_S= (L_SS+mu*(M_SS-eye(n_S)/2))*sigma;
   v_S=(Mt_SS+eye(n_S)/2+mu*N_SS)*sigma;
  
% calculate L_DS, M_DS, the discrete form of theoperators for the domain points
%   dummy values set to vecp_D, since it is not used
   for (i=1:n_D)
     vecp_D(i,1)=1;
     vecp_D(i,2)=0;
   end   
   [L_DS,M_DS,Mt_DS,N_DS] =lbem2(n_D,p_D,vecp_D,n_S,vertpts_S,elemvert_S,false,true,true,false,false);
   
   
% calculate the solution at the domain points
   phi_D=(L_DS+mu*M_DS)*sigma;