%***************************************************************
% interiorSquareTestRobin.m by www.boundary-element-method.com *
% **************************************************************
% Matlab/Freemat/Octave codes
% interiorSquareTestRobin.m Version 1 July 2015
% Stephen Kirkup
% School of Engineering, University of Central Lancashire, UK
% http://www.researchgate.net/profile/Stephen_Kirkup
% Issued under the GNU General Public License 2007
%
% Source of the code: 
%  http://www.boundary-element-method.com/mfiles/
%   interiorSquareTestRobin.m
% Source of the user-guide: 
%  http://www.boundary-element-method.com/mfiles/
%   interiorSquareTestRobin_m.pdf
% Code tested: 
%  http://www.boundary-element-method.com/mfiles/gls.m
%  http://www.boundary-element-method.com/mfiles/regls.m

% www.boundary-element-method.com
% contact: stephen@boundary-element-method.com

% function [phi_D,phi_S,v_S] = interiorSquareTestRobin
% Returns the solution phi_D, the solution at the domain points, phi_S and v_S,
% the solution on the boundary. The boundary is defined internally by
% calling 'square' and the points are defined by calling 'points'.

function [phi_D,phi_S,v_S]= interiorSquareTestRobin(n_S)
    
    format long
    if (rem(n_S,4)>0)
        'n_S must be a multiple of 4. Program halted'
        halt
    end
        if (n_S<4)
        'n_S must be 4 or greater'
        halt
    end
        
    
% Set the square surface
%  The length of the sides of the square
    dimension=1.0;

% Declare the following as column vectorsby setting each as a column of zeros 
    v_S=zeros(n_S,1);
    phi_S=zeros(n_S,1);
    c=zeros(n_S,1);
    alpha=zeros(n_S,1);
    beta=zeros(n_S,1);
    f=zeros(n_S,1);
    f_new=zeros(n_S,1);
    phi_S_exact=zeros(n_S,1);
    v_S_exact=zeros(n_S,1);
    
%  Create the mesh
    centre=[0.5 0.5];
    dimension=1;
    clockwise=1;
    [vertpts,elemvert]=square_general(n_S,dimension,centre,clockwise);
%  Set the points either side of each panel 
    [qa,qb]=vertices(n_S,vertpts,elemvert);
%  Set the collocation points, at the centre of each panel
    colpoints=(qa+qb)/2;
%  Set the 5 interior points
    n_D=2;
    p_D=[ 0.25, 0.75;  0.75, 0.25];
   
%  Set the matrices corresponding to the Laplace Integral Operators for the
%   collocation points on the boundary  
    [L_SS,M_SS,Mt_SS,N_SS] =lbem2_on(n_S,vertpts,elemvert,true,true,true,false);
   
%  Set the matrices corresponding to the Laplace Integral Operators for the
%   points in the domain
%   Dummy values are set to vecp_D, since it is not used
    for (i=1:n_D)
      vecp_D(1,i)=1;
      vecp_D(2,i)=0;
    end
%   Calculate L_DS, M_DS, Mt_DS and N_DS, the discrete form of the operators for the domain points
    [L_DS,M_DS,Mt_DS,N_DS] =lbem2(n_D,p_D,vecp_D,n_S,vertpts,elemvert,false,true,true,false,false);

% The M_SS and Mt_SS matrices are only used in these fomulations with I/2 added
    M_SSplus=M_SS+eye(n_S)/2;
    Mt_SSplus=Mt_SS+eye(n_S)/2;

    
% State the number of elements    
   'number of elements'
   n_S
   ' '
% For this set of tests {\phi}=x^2+y^2
%  Set the Robin boundary condition 
    for k=1:n_S
        phi_S_exact(k)=2*(colpoints(k,1)^2-colpoints(k,2)^2);
    end
    for k=1:n_S/4
        v_S_exact(k)=0;
        v_S_exact(n_S/4+k)=-4;
        v_S_exact(n_S/2+k)=4;
        v_S_exact(3*n_S/4+k)=0;
    end
    
% For this set of tests {\phi}=x^2+y^2
%  Set the Robin boundary condition - half Dirichlet and half Neumann
    alpha(1:n_S/2)=1.0;
    alpha(n_S/2+1:n_S)=0.0;
    beta(1:n_S/2)=0.0;
    beta(n_S/2+1:n_S)=1.0;
    for (k=1:n_S)
       f(k)=alpha(k)*phi_S_exact(k)+beta(k)*v_S_exact(k);
    end
 
 
    '1 Direct Robin L v = (M+I/2){\phi}'
%Solution using row swaps
% Set the matrices so that A*sol=B*F'
    A=[ L_SS(1:n_S/2,1:n_S/2)                    -M_SS(1:n_S/2,n_S/2+1:n_S);
    L_SS(n_S/2+1:n_S,1:n_S/2)                -M_SSplus(n_S/2+1:n_S,n_S/2+1:n_S)];
    B=[ M_SSplus(1:n_S/2,1:n_S/2)   -L_SS(1:n_S/2,n_S/2+1:n_S);
         M_SS(n_S/2+1:n_S,1:n_S/2)               -L_SS(n_S/2+1:n_S,n_S/2+1:n_S)];
% Solve the system        
    sol=A\(B*f);
% Complete v_S and phi_S
    v_S(1:n_S/2)=sol(1:n_S/2);
    v_S(n_S/2+1:n_S)=f(n_S/2+1:n_S);
    phi_S(1:n_S/2)=f(1:n_S/2);
    phi_S(n_S/2+1:n_S)=sol(n_S/2+1:n_S);
% Find phi_D
    phi_D=L_DS*v_S-M_DS*phi_S

% Solution using gls    
    [phi_S, v_S, xory, L, U, B_gls, P, lfail ] = gls(M_SSplus, L_SS, c, n_S, alpha, beta, f);    
% Find phi_D
    phi_D=L_DS*v_S-M_DS*phi_S
    
% Solution using 2nx2n matrix
    phi_v_S= [M_SSplus -L_SS; diag(alpha) diag(beta)]\[c;f];
    phi_S=phi_v_S(1:n_S);
    v_S=phi_v_S(n_S+1:2*n_S);
% Find phi_D
    phi_D=L_DS*v_S-M_DS*phi_S

    '1 Indirect Robin {\phi}=L{\sigma} on S and in D, v=(Mt-I/2){\sigma}'
% Set the matrices so that A*sigma=f
    for i=1:n_S
        for j=1:n_S
            A(i,j)=alpha(i)*L_SS(i,j)+beta(i)*Mt_SSplus(i,j);
        end
    end
    
    sigma_S=A\f;
% Complete v_S and phi_S
    phi_D= L_DS*sigma_S

% For this set of tests {\phi}=x^2+y^2
%  Set the Robin boundary condition 
    for k=1:n_S
        alpha(k)=10^(6*k/n_S);
        beta(k)=10^(6*(n_S-k)/n_S);
        f(k)=alpha(k)*phi_S_exact(k)+beta(k)*v_S_exact(k);
    end        

    '2 Direct Robin L v = (M+I/2){\phi}'    
% Solution using gls    
    [phi_S, v_S, xory, L, U, B_gls, P, lfail ] = gls(M_SSplus, L_SS, c, n_S, alpha, beta, f);    
% Find phi_D
    phi_D=L_DS*v_S-M_DS*phi_S
    for k=1:n_S
        f_new(k)=2*f(k);
    end         
    [phi_S2,v_S2]=regls(n_S, xory, L, U, B_gls, P, alpha, beta, c, f_new);  
    phi_D2=L_DS*v_S2-M_DS*phi_S2

  
    
% Solution using 2nx2n matrix
    phi_v_S= [(M_SS+eye(n_S)/2.0) -L_SS; diag(alpha) diag(beta)]\[c;f];
    phi_S=phi_v_S(1:n_S);
    v_S=phi_v_S(n_S+1:2*n_S);
% Find phi_D
    phi_D=L_DS*v_S-M_DS*phi_S

    '2 Indirect Robin {\phi}=L{\sigma} on S and in D, v=(Mt-I/2){\sigma}'
% Set the matrices so that C*sigma=f
    for i=1:n_S
        C(i,1:n_S)=alpha(i)*L_SS(i,1:n_S)+beta(i)*Mt_SSplus(i,1:n_S);
    end   
    sigma_S=C\f;
% Complete v_S and phi_S
    phi_D= L_DS*sigma_S


     
end
      
      
      
      
    
  
     
  
  
           
