
function demo_nurbs_cubature_02(example)

%---
% OBJECT:
%----
% Demo illustrating CUBATURE:
% 1. how to define a NURBS on a composite boundary, using arcs of disks,
%    ellipses, segments, polygons and "free NURBS". The routines work on
%    piecewise NURBS of different order.
%    To this purpose see the MATLAB function "define_domain" at the
%    bottom of this demo.
% 2. how to get an algebraic cubature rule on the desired domain, via
%    "cubRS" routine.
% 3. testing the routines over random polynomials.
%---
%% Copyright (C) 2007-2024 Chiara Zara, Alvise Sommariva, Marco Vianello.
%%
%% 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 2 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, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Author:  
%%          Chiara Zara
%%          Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: July 8, 2024 -
%---
% Dates
%---
% First version: October 31, 2021;
% Checked: July 8, 2024.
%---


warning off;

% ------------------------------ Settings ---

% Algebraic degree of precision of the cubature rule. Suggested values are
% mild, e.g. not larger than "ade=15". Typical "ade" in the applications
% are between "3" and "5".
adeV=2:2:10;

% Number of tests to determine the median cputime. We suggest to set
% "tests=10", while we do not recommend "tests=1" because it often turns
% out to provide unreliable cputimes.
tests=1;

% Cubature extraction method: the variable "extraction_type" sets how
% Tchakaloff rule is extracted.
%      1: lsqnonneg, 2: lawsonhanson, 3: LHDM 4: \
% If not declared or equal to "[]", the default is LHDM, see
% "dCATCH: a numerical package for d-variate near G-optimal Tchakaloff
% regression via fast NNLS"
% M. Dessole, F. Marcuzzi and M. Vianello
% MDPI-Mathematics 8 (2020) - Special Issue "Numerical Methods"
extraction_type=3;

% Nbox is a technical parameter for indomain routine; in doubt set 100.

Nbox=3;

% In indomain routine one must be certain that points are inside the
% domain. In our cubature needs we just wants some points in the domain.
% Consequently, if we are not fully sure that a point is in the domain, due
% to geometrical issues, we do not take it into account.

safe_mode=1;

% The variable "domain_type" set the domains taken into account.
if nargin < 1
    example=1;
end









% ----------------------------- Main code ---

%---
% 1. Make NURBS structure
%---
[geometry_NURBS_EXT,geometry_NURBS_INT]=define_domain(example);

REV=[]; logerrV=[];

switch example
    case 1
        x0=0; y0=0.4;
    case {2,4}
        x0=-0.2; y0=-1.4;
    case 3
        x0=0.5; y0=-1;
end

% Define functions
f1=@(x,y) exp(-(x.^2+y.^2));
f2=@(x,y) ((x-x0).^2+(y-y0).^2).^(11/2);
f3=@(x,y) ((x-x0).^2+(y-y0).^2).^(1/2);

% Define high order rule
ADEref=30;
xywL = cubRS(ADEref,geometry_NURBS_EXT,geometry_NURBS_INT,0,Nbox,safe_mode);

xxL=xywL(:,1); yyL=xywL(:,2); wwL=xywL(:,3);
I(1)=wwL'*feval(f1,xxL,yyL);
I(2)=wwL'*feval(f2,xxL,yyL);
I(3)=wwL'*feval(f3,xxL,yyL);

Zcards=[];

for jj=1:length(adeV)

    ade=adeV(jj);

    %---
    % 2. Compute algebraic cubature rule of the domain (several tests!).
    %---
    [xyw,res,Z,Zin] = cubRS(ade,geometry_NURBS_EXT,geometry_NURBS_INT,...
        extraction_type,Nbox,safe_mode);

    Zcards(end+1,:)=[size(xyw,1) size(Zin,1) size(Z,1)];

    xx=xyw(:,1); yy=xyw(:,2); ww=xyw(:,3);

    %---
    % 3. Making tests.
    %---
    I0(1)=ww'*feval(f1,xx,yy);
    I0(2)=ww'*feval(f2,xx,yy);
    I0(3)=ww'*feval(f3,xx,yy);

    RE=abs((I-I0)./I);

    %---
    % 4. Statistics.
    %---

    fprintf('\n \t %2.0f | %1.1e | %1.1e | %1.1e | ',ade,RE(1),RE(2),RE(3));


end

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

for k=1:size(Zcards,1)
    fprintf('\n \t %2.0d | %4d | %4d  | %4d  | ',adeV(k),Zcards(k,1),...
        Zcards(k,2),Zcards(k,3)); 
end

fprintf('\n \n');








%---
% 6. Plots.
%---

h=figure(1);
f1=ishandle(h)&&strcmp(get(h,'type'),'figure'); if f1,clf(1);end
axis equal;
axis tight;
plotNURBSPL(geometry_NURBS_EXT);
hold on;
plotNURBSPL(geometry_NURBS_INT);
plot(xx,yy,'mo','MarkerEdgeColor','k',...
                       'MarkerFaceColor','m',...
                       'MarkerSize',4)      
hold off;























