solverEllipsoidPenetrable.m

% Solve acoustic scattering for a penetrable ellipsoid
%
%  s = solverEllipsoidPenetrable(k,inc,x0,ar,refin,rho,n) initialises a solver class to
%  simulate acoustic scattering by an ellipsoid using the method of
%  fundamental solutions. Parameters are:
%
%    k - the wavenumber
%    inc - specifies the incident wave (of class incident or empty)
%    x0 - the centre of the ellipsoid (a 3 x 1 vector)
%    ar - the radius in the x,y,z directions (a 3 x 1 vector)
%    refin - the relative refractive index of the ellipsoid
%    rho - a vector [rho0 rho1] where rho0 and rho1 are the density outside
%    and inside the ellipsoid respectively
%    n - parameter controlling the number of discrete sources in the interior
%
%  s = solverEllipsoidSoft(k,inc,x0,ar,refin,rho,n,m) further specifies:
%
%    m - parameter controlling the number of surface points at which the
%    boundary condition is matched.
%
% Also:
%
%  s.setup() prepares the class. This only needs to be done once and is
%  independent of the incident wave.
%
%  s.solve() solves the scattering problem. The object s can be
%  subsequently used to compute the far field, cross section etc.
%
%  s.visualise() visualises the scatterer.
%
%  s.setIncidentField(inc) sets the incident field as specified in the cell
%  array inc.
%
%  val = s.getFarField(pts) computes the far field for the
%  incident field inc{1} at points on the unit sphere specified in a 3 x n 
%  matrix pts.
%
%  val = s.getFarField(pts,index) computes the far field for the
%  incident fields inc{index} at points on the unit sphere specified in a 
%  3 x n matrix pts.
%
%  val = s.getCrossSection(pts) computes the cross section for the
%  incident field inc{1} at points on the unit sphere specified in a 3 x n 
%  matrix pts.
%
%  val = s.getCrossSection(pts,index) computes the cross section for the
%  incident fields inc{index} at points on the unit sphere specified in a 
%  3 x n matrix pts.
%  
% See also: solver, solverEllipsoidSoft, solverEllipsoidHard.
%
% Stuart C. Hawkins - 21 April 2021

% Copyright 2019-2022 Stuart C. Hawkins
% 	
% This file is part of TMATROM3
% 
% TMATROM3 is free software: you can redistribute it and/or modify	
% it under the terms of the GNU General Public License as published by	
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% TMATROM3 is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with TMATROM3.  If not, see <http://www.gnu.org/licenses/>.


classdef solverEllipsoidPenetrable < solver

    properties
        origin
        shell
        basis_minus
        basis_plus
        matrix
        boundary_conditions
        cof
        basis_param
        ar
        refin
        rho
    end
    
    methods
        

        %-----------------------------------------
        % constructor
        %-----------------------------------------
        
        function self = solverEllipsoidPenetrable(kwave,incidentField,origin,ar,refin,rho,basis_param,m)

            % Note: rho = [rho_0 rho_1]
            
            % set default for m
            if nargin<8
                m = basis_param;
            end
            
            % call parent constructor
            self = self@solver(kwave,incidentField);
            
            % store aspect ratio
            self.ar = ar;
            
            % store origin
            self.origin = origin;
            
            % store refractive index and density
            self.refin = refin;
            self.rho = rho;

            % store number of sources parameter
            self.basis_param = basis_param;

            % setup boundary
            bdry = translated_sheet(self.origin(:),ellipsoid_segment([0 pi],[0 2*pi],ar(1),ar(2),ar(3)));
            self.shell = sheet_with_points(bdry,midpoint_rule(m),midpoint_rule(2*m));
            
            % basis for exterior
            int = translated_sheet(self.origin(:),ellipsoid_segment([0 pi],[0 2*pi],0.95*ar(1),0.95*ar(2),0.95*ar(3)));
            minus = sheet_with_points(int,midpoint_rule(basis_param),rectangle_left_rule(2*basis_param));
            self.basis_minus = basis(minus,kwave);
            
            % basis for interior
            ext = translated_sheet(self.origin(:),ellipsoid_segment([0 pi],[0 2*pi],1.05*ar(1),1.05*ar(2),1.05*ar(3)));
            plus = sheet_with_points(ext,midpoint_rule(basis_param),rectangle_left_rule(2*basis_param));
            self.basis_plus = basis(plus,refin*kwave);

            % record boundary conditions
            self.boundary_conditions = ...
                [dirichlet_relation([-1 1],[self.basis_minus self.basis_plus],self.shell,[]);...
                neumann_relation([-1 rho(1)/rho(2)],[self.basis_minus self.basis_plus],self.shell,[])];
            
        end
                 
        %-----------------------------------------
        % setup
        %-----------------------------------------
        
        function setup(self)                      
            
            % compute the system matrix
            self.matrix = self.boundary_conditions.matrix();
            
        end
        
        %-----------------------------------------
        % solve
        %-----------------------------------------
        
        function solve(self)

            % loop through the right hand sides to assemble the right hand
            % side matrix
            for j=1:length(self.incidentField)
            
                % associate jth incident field with the boundary conditions
                % Note: for penetrable there are two boundary conditions...
                % (1) is the Dirichlet and (2) is the Neumann
                self.boundary_conditions(1).fun = self.incidentField{j};
                self.boundary_conditions(2).fun = self.incidentField{j};
            
                % evaluate the vector associated with the incident field
                rhs(:,j) = self.boundary_conditions.rhs();
                
            end
            
            % Matlab sometimes warns about a rank deficient matrix but this
            % can be safely ignored so turn it off
            warning('off','MATLAB:rankDeficientMatrix')

            % solve the linear system
            % Note: self.matrix is assembled in the setup() method
            self.cof = self.matrix \ rhs;
            
            % turn the warning back on
            warning('on','MATLAB:rankDeficientMatrix')
            
        end
        
        %-----------------------------------------
        % get far field
        %-----------------------------------------
        
        function val = getFarField(self,points,index)

            % set default for index
            if nargin<3
                index = 1:length(self.incidentField);
            end
            
            % loop through index
            for j=1:length(index)
                
                % compute the far field associated with right hand side
                % index(j)
                val(:,j) = self.basis_minus.farfield(points,self.cof(1:2*self.basis_param^2,j));
                
            end            

        end
        
        %-----------------------------------------
        % text description of the solver
        %-----------------------------------------
        
        function str = description(self)
           
            str = sprintf('Method of Fundamental Solution solver for penetrable ellipsoid (refractive index %0.1f and density %0.1f) with semi-axes (%0.1f,%0.1f,%0.1f).',...
                self.refin,self.rho(2),self.ar(1),self.ar(2),self.ar(3));
            
        end
        
        %===============================================================
        % you may provide other methods required to implement your solver
        % or help the user
        %===============================================================

        %-----------------------------------------
        % visualise the scatterer
        %-----------------------------------------

        function varargout = visualise(self)
        
            % number of points in each dimension in the surface mesh
            n = 60;
            
            % setup mesh points in the reference rectangle
            u = linspace(0,1,n);
            v = linspace(0,1,n);           
            [uu,vv] = meshgrid(u,v);

            % transform the mesh points to the surface
            % Note: sheet.evaluate() method wants points in vectors so we
            % turn uu and vv into column vectors first
            f = self.shell.sheet.evaluate(uu(:),vv(:));
            
            % reshape the result into an n x n matric, corresponding to the
            % grid
            f = reshape(f,3,n,n);
                  
            % surf the surface
            shell = surf(squeeze(f(1,:,:)),squeeze(f(2,:,:)),squeeze(f(3,:,:)));
            
            % tweak the appearance
            set(shell,'facecolor',0.5*[1 1 1])
            set(shell,'edgecolor','none')            

            % return handle of surface if required
            if nargout > 0
                varargout{1} = shell;
            end
            
        end
        
    end % end methods
    
end