
function demo_nurbs_cubature_01(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=100;

% 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=100;

% 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 domain taken into account. 
if nargin < 1
  example=3;
end









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

%---
% 1. Make NURBS structure
%----

[geometry_NURBS,geometry_NURBSH,domain_str]=define_domain(example);

REV=[]; logerrV=[];
for jj=1:length(adeV)

    ade=adeV(jj);

    %---
    % 2. Compute algebraic cubature rule of the domain (several tests!).
    %---
    xyw = cubRS(ade,geometry_NURBS,geometry_NURBSH,extraction_type,Nbox,...
        safe_mode);

    %---
    % 3. Compute algebraic cubature rule of the boundary.
    %---

    [inside,in_doubts,sp_xnV,sp_xdV,sp_ynV,sp_ydV,boxVx,...
        singular_points_X,singular_points_Y]=inRS([],geometry_NURBS,...
        Nbox,safe_mode);
    xywL = lineint(ade+1,sp_xnV,sp_xdV,sp_ynV,sp_ydV);

    [insideH,in_doubtsH,sp_xnVH,sp_xdVH,sp_ynVH,sp_ydVH,boxVxH,...
        singular_points_XH,singular_points_YH]=inRS([],geometry_NURBSH,...
        Nbox,safe_mode);
    xywLH = lineint(ade+1,sp_xnVH,sp_xdVH,sp_ynVH,sp_ydVH);

    %---
    % 4. Making tests.
    %---

    for k=1:tests
        a=rand(1); b=rand(1); c=rand(1);
        f=@(x,y) (a+b*x+c*y).^ade;

        xx=xyw(:,1); yy=xyw(:,2); ww=xyw(:,3);
        I0=ww'*feval(f,xx,yy);

        xxL=xywL(:,1); yyL=xywL(:,2); wwL=xywL(:,3);
        g=@(x,y) (a+b*x+c*y).^(ade+1)/(b*(ade+1));
        I1E=wwL'*feval(g,xxL,yyL);

        xxLH=xywLH(:,1); yyLH=xywLH(:,2); wwLH=xywLH(:,3);
        I1H=wwLH'*feval(g,xxLH,yyLH);

        I1=I1E-I1H;

        RE(k,1)=abs((I0-I1)./I1);
    end

    %----
    % 5. Statistics.
    %----
    fprintf('\n \n \t ADE   : %2.0f',ade);
    fprintf('\n \t RE max: %1.3e',max(RE));
    kpos=find(RE > 0);
    logerr=10^(sum(log10(RE(kpos)))/length(kpos));
    fprintf('\n \t RE log: %2.3e',logerr);

    cell_tick{jj}=num2str(ade);

    REV=[REV RE];
    logerrV=[logerrV; logerr];
end

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

h=figure(1);
f1=ishandle(h)&&strcmp(get(h,'type'),'figure'); if f1,clf(1);end
figure(1)
axis equal;
for jj=1:length(adeV)
    n=adeV(jj);
    reL=REV(:,jj); log_reL=logerrV(jj);
    plot_errors(jj,n,reL,log_reL);
    hold on;
end
ax=gca;
ax.XAxis.FontSize = 16;
ax.YAxis.FontSize = 16;
xlim([adeV(1)-1,adeV(end)+1]);
xticks(adeV)
hold off
fprintf('\n \n');

h=figure(2);
f2=ishandle(h)&&strcmp(get(h,'type'),'figure'); if f2,clf(2);end
figure(2)

axis equal;
plotNURBSPL(geometry_NURBS);
hold on;
plotNURBSPL(geometry_NURBSH);
plot(xx,yy,'mo','MarkerEdgeColor','k',...
                       'MarkerFaceColor','m',...
                       'MarkerSize',4)      
hold off;
















