
function demo_glaubitz_S2

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Demo used to compute integrals from samples involving scattered data.
%
% 1. Determines domain and functions to study from function defined by the
%    user.
% 2. computes integral using only evaluations at scattered data:
%
% (a) For methods based on evaluating in scattered data, using the values
%    to evaluate at points of an algebraic rule:
%    Makes experiments, evaluating the integrand at cubature nodes of an
%    algebraic rule, when only samples at scattered data are available.
%
%    Important: algebraic degree of exactness of the cubature rule is fixed.
%
% (b) By Glaubitz algorithm computes the value of the integral, with a 
%    formula based on scattered data.
%
% 3. Makes plot of the domain, scattered samples, cubature nodes.
%--------------------------------------------------------------------------
% Modified: November 6, 2024.
%--------------------------------------------------------------------------
% COPYRIGHT
%--------------------------------------------------------------------------
% Copyright (C) 2024-
%
% Authors:
% Alvise Sommariva.
%
% This program 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.
%
% This program 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 this program.  If not, see <https://www.gnu.org/licenses/>.
%--------------------------------------------------------------------------



% ...................... Problem settings  ..........................

%--------------------------------------------------------------------------
% 
%--------------------------------------------------------------------------
% Note: use deg_rule <=15 for NURBS
%--------------------------------------------------------------------------

% 1: 'sphere';
% 2: 'spherical-triangle';
% 3: 'spherical-rectangle';
% 4: 'spherical-polygon';
domain_type=4;

card=800;    % scattered data cardinality
scat_type='halton'; % scattered data type.

f=@(x,y,z) 1./( (1+x.^2).*(1+y.^2).*(1+z.^2));

method='Glaubitz-algorithm';




% ....................... Preliminary settings  ...........................

% some additional routines are required
addpath(genpath('../EXTERNAL_ROUTINES/'));

% flag_compression: if "1" compression will be attempted, if necessary,
% otherwise the code will propose an alternative rule.
flag_compression=1;



% ...................... Making numerical tests  ..........................


% 1. Define domain
%    Note: the domain routine is in "../EXTERNAL_ROUTINES/DCUB" folder

fprintf('\n \t * defining integration domain');
domain_example=domain_str(domain_type);
domain_struct=define_domain(domain_example);
fprintf('\n \n \t \t The domain is: '); disp(domain_example);

% 1A. Define scattered data
fprintf('\n \t * defining scattered data');
[t,dbox,area_domain]=define_scattered_pointset(card,domain_struct,...
    scat_type);
ft=feval(f,t(:,1),t(:,2));

% 1B. Define cubature rule
%    Note: the domain routine is in "../EXTERNAL_ROUTINES/DCUB" folder
XYW=define_cub_rule(domain_struct,deg_rule,flag_compression);

% 1C: Compute integral from evaluation at scattered data


area_dbox=diff(dbox(1,:))*diff(dbox(2,:));
wQMC=(area_domain/area_dbox)/size(t,1);
[w,deg_rule]=glaubitz_algorithm(t,domain_struct,wQMC,'LS');
I=w'*ft;


%--------------------------------------------------------------------------
% 2. Display statistics.
%--------------------------------------------------------------------------
fprintf('\n \t -------------------------------------------------------');
fprintf('\n \t Domain (integration)          : '); disp(domain_struct.domain);
fprintf('\n \t Scattered data type           : '); disp(scat_type);
fprintf(   '\t Evaluation method             : '); disp(method);
fprintf('\n \t Scatt. points card. required  : %5.0f',card);
fprintf('\n \t Scatt. points card. provided  : %5.0f \n',size(t,1));

fprintf('\n \t Degree of precision           : %5.0f',deg_rule);

card_rule=size(XYW,1);
fprintf('\n \t Cardinality of the formula    : %5.0f',card_rule);

dim_poly=(deg_rule+1)*(deg_rule+2)/2;
fprintf('\n \t Dimension of the poly. space  : %5.0f',dim_poly);


fprintf('\n \t Integral approximation        : %1.15e',I);

switch domain_type
    case 17
        atol=10^(-12); rtol=10^(-12);
        IR = integral2(f,-1,1,-1,1,'AbsTol',atol,'RelTol',rtol);
    case 18
        atol=10^(-12); rtol=10^(-12);
        IR = integral2(f,0,1,0,1,'AbsTol',atol,'RelTol',rtol);
    otherwise
        ft=feval(f,XYW(:,1),XYW(:,2));
        IR=(XYW(:,3))'*ft;
end

fprintf('\n \t Integral computed by alg.rule : %1.15e',IR);
fprintf('\n \t Absolute error                : %1.1e',abs(IR-I));
fprintf('\n \t Relative error                : %1.1e',abs(IR-I)/abs(I));

fprintf('\n \n');

fprintf('\n \t ........................................................ ');

fstr = strrep(char(f),'@(x)','');
fprintf('\n \n \t The integrand is: '); disp(fstr);
fprintf('\n \t The domain is   : '); disp(domain_struct.domain);
fprintf('\n \n \t ........................................................ ');
fprintf('\n \n');



% 3. plot domain and scattered data
figure(1)
plot_2D(domain_struct,t,XYW(:,1:2));
saveas(gcf,domain_example,'epsc');









function domain_example=domain_str(example)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Given a number "example", it provides the name of the corrsponding domain
% in the gallery,
%--------------------------------------------------------------------------

switch example
    case 1, domain_example='sphere';
    case 2, domain_example='spherical-triangle';
    case 3, domain_example='spherical-rectangle';
    case 4, domain_example='spherical-polygon';
end

















