% BEM-LAP-MAT project
% Matlab/Freemat codes
% Copyright 2008 Stephen Kirkup 
% John Tyndall Research Institute
% UCLan Westlakes Campus, Samuel Lindow Building, Westlakes Science and Technology Park, Whitehaven, CA24 3JY
% Issued under the GNU General Public License 2007, see gpl.txt
% www.boundary-element-method.com
% www.kirkup.info/opensource


% function [phi_D,phi_S,v_S] = libem2(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(n_D,p_D,n_S,vertpts_S,elemvert_S,alpha_S,beta_S,f_S)
   
  % calculate phi_S, v_S, phi and v on the boundary

  %  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);
 
   M_SSplus=M_SS+eye(n_S)/2;
   Mt_SSminus=Mt_SS-eye(n_S)/2;
   mu=norm(L_SS)/norm(Mt_SSminus);
   for (i=1:n_S)
      zero_S(i)=0.0;
   end
   [phi_S,v_S]=gls(M_SSplus+mu*N_SS,L_SS+mu*Mt_SSminus,zero_S,n_S,alpha_S,beta_S,f_S);
   
  % calculate phi_S, v_S, phi and v on the boundary
   
  %  calculate L_DS, M_DS, Mt_DS and N_DS, the discrete form of the operators 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);
  
  phi_D=L_DS*v_S-M_DS*phi_S;