
function [XYZW,flag,iteration_stats]=cubature_polyhedron_free(...,
    vertices,facets,ade,extraction_type,set_type)

%--------------------------------------------------------------------------
% Object:
%
% This routine computes nodes and weights of a cubature formula with degree
% of precision "ade" on a polyhedron defined by "vertices" and "facets".
%
% The cardinality of the set stored in XYZW is at most 
%             "(ade+1)*(ade+2)*(ade+3)/6" 
% and the nodes are inside the domain.
%
% A first attempt is made to find a rule with positive weights, otherwise
% it is found one with some possible negative weights, but in general
% "stable".
%--------------------------------------------------------------------------
% Input:
%
% vertices: The vertices are stored in "vertices" and is an "m x 3" matrix.
%  The k-th vertex has coordinates "vertices(k,:)".
%
% facets: The variable "facets" can be:
%
% 1. A "cell array":
%    In this case the i-th cell is a vector of positive integers,
%    representing the indices of vertices of the i-th facet, taking into
%    account the variable "vertices". The orientation of the facet is
%    counterclockwise (w.r.t to outward normal).
%
% 2. A "matrix" with dimension "m x 3":
%    In this case the i-th row is a vector of positive integers,
%    representing the indices of vertices of the triangular facet, taking
%    into account the variable "vertices". The orientation of the
%    triangular facet is counterclockwise (w.r.t to outward normal).
%
% ade: algebraic degree of precision of the rule (total degree).
%
% * extraction_type: algorithm that implements compression:
%    0: Dessole, Marcuzzi, Vianello algorithm (default)
%    1. lsqnonneg
%    2. Slawski's Lawson-Hanson algorithm
%    3. QR compression
%
% * set_type: the nodes are taken from certain pointsets:
%    1. tensorial type
%    2. Halton type (default)
%    3. random type
%
% Note: The variables in which there is an asterisk are not mandatory.
%--------------------------------------------------------------------------
% Output:
%
% XYZW: The variable "XYZW" is a matrix with "m" rows (that is the
%  cardinality of the cubature pointset) and 4 columns.
%  The first 3 columns are the coordinates of the nodes of a
%  rule over a reference tetrahedron, while the 4th column represent the
%  corresponding weights.
%
% flag: 0: termination with failure, computing rule with possible negative
%          weights and internal nodes
%       1: termination with success, computing rule with positive weights
%          and internal nodes.
%
% iteration_stats: the method perform several iteration and in each one a
%      pointset is analysed so to find the wanted cubature rule that
%      hopefully matches the shifted-tensorial chebyshev moments of total
%      degree "ade"; this variable is a L x 4 matrix in which each rows
%      describes an iteration and is of the form
%
%             [iteration full_set_card inner_set_card residual iters]
%
%     where:
%
%     * full_set_card: cardinality of the pointset generated in an
%                      hyperrectangle containing the domain;
%
%     * inner_set_card: cardinality of the points of the previous pointset
%                      that are in the polyhedron;
%
%     * residual: how the generated cubature provides moments that are
%                 close to those of the shifted-tensorial chebyshev basis,
%                 based on the infinity norm of their difference.
%
%     * iters: in the case of Dessole-Marcuzzi-Vianello algorithm or Matlab
%     built-in lsqnonneg it gives the number of inner iterations of the
%     method.
%
%--------------------------------------------------------------------------
% Routines called by "cubature_polyhedron_tenscheb_02".
%
% 1. determine_facets_planes (attached below)
% 2. test_convexity (attached below)
% 3. chebyshev_moments_polyhedron (and its subroutines, external)
% 4. monomials_moments_convex_polyhedron (and its subroutines, external)
% 5. rule_search_chebyshev_01 (and its subroutines, main routine is
%    attached below).
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------

% ........................... TROUBLESHOOTING  ............................

if nargin < 4, extraction_type=0; end
if nargin < 5, set_type=2; end

% ........................... CONVEXITY CHECK .............................

[A,b]=determine_facets_planes(vertices,facets);
isconvex=test_convexity(A,b,vertices);


% ........................... COMPUTE MOMENTS .............................

moments_ch=chebyshev_moments_polyhedron(vertices,facets,ade);
vol_poly=moments_ch(1); % polyhedron volume

% .............................. RULE SEARCH ..............................

[XYZW,flag,iteration_stats]=rule_search_chebyshev_01(ade,vertices,...
    facets,moments_ch,extraction_type,set_type,isconvex,A,b,vol_poly);









function flag=test_convexity(A,b,points)

%--------------------------------------------------------------------------
% Object:
% Testing polyhedron convexity.
%--------------------------------------------------------------------------
% Input:
%
% A: matrix N x 3 describing the facets hyperplanes so that if the domain
% is convex then it satisfies A*[X Y Z]' <= b for each of its point [X Y Z]
%
% b: vector N x 1, determining the domain (if convex)
%
% points: matrix M x 3 (polyhedron vertices coordinates as row vector)
%
%--------------------------------------------------------------------------
% Output:
% flag: 0: not convex 1: convex
%--------------------------------------------------------------------------

c=A*points';
tol=10^(-14);
points_test= (c <= b+tol);

L=size(A,1);
sum_points_test=norm(sum(points_test,1)-L);
if sum_points_test == 0
    flag=1;
else
    flag=0;
end









function [XYZW,flag,iteration_stats]=rule_search_chebyshev_01(n,vertices,...
    facets,moments,extraction_type,set_type,isconvex,A,b,volume_poly)


%--------------------------------------------------------------------------
% Object:
%
% This routine extracts nodes and weights of a cubature formula with degree
% of precision "ade" on a polyhedron defined by "vertices" and "facets".
% The cardinality of the set is at most "(ade+1)*(ade+2)*(ade+3)/6" and the
% nodes are inside the domain.
% A first attempts is made to find a rule with positive weights, otherwise
% it is found one with some possible negative weights, but in general
% "stable".
%
%--------------------------------------------------------------------------
% Input:
% vertices: The vertices are stored in "vertices" and is an "m x 3" matrix.
%  The k-th vertex has coordinates "vertices(k,:)".
%
% facets: The variable "facets" can be:
%
% 1. A "cell array":
%    In this case the i-th cell is a vector of positive integers,
%    representing the indices of vertices of the i-th facet, taking into
%    account the variable "vertices". The orientation of the facet is
%    counterclockwise (w.r.t to outward normal).
%
% 2. A "matrix" with dimension "m x 3":
%    In this case the i-th row is a vector of positive integers,
%    representing the indices of vertices of the triangular facet, taking
%    into account the variable "vertices". The orientation of the
%    triangular facet is counterclockwise (w.r.t to outward normal).
%
% moments: moments of a shifted-tensorial chebyshev basis on polyhedron,
%    say "T_i(x)*T_j(y)*T_k(z)".
%    The shift is based on the smaller hyperrectangle containing the
%    domain, having sides parallel to the x-y-z axis.
%    The moments are ordered lexicographically w.r.t. (i,j,k).
%
% ade: algebraic degree of precision of the rule (total degree).
%
% extraction_type: algorithm that implements compression:
%    0: Dessole, Marcuzzi, Vianello algorithm (default)
%    1. lsqnonneg
%    2. Slawski's Lawson-Hanson algorithm
%    3. QR compression
%
% set_type: the nodes are taken from certain pointsets:
%    1. tensorial type
%    2. Halton type (default)
%    3. random type
%
% isconvex: 0 if the domain is convex, 1: otherwise.
%   if the polyhedron D is convex, the needed indomain routine is based on
%   the fact that X belongs to D if and only if A*X <= b (for certain A,b).
%
% A, b: matrix and vector such that if the polyhedron D is convex, then a
%    point X belongs to D if and only if A*X <= b (for certain A,b).
%
% volume_poly: approximation of the volume of the domain (it may be unknown
%    and in which case do not set this variable or make it equal to NaN).
%
%--------------------------------------------------------------------------
% Output:
%
% XYZW: The variable "XYZW" is a matrix with "m" rows (that is the
%  cardinality of the cubature pointset) and 4 columns.
%  The first 3 columns are the coordinates of the nodes of a
%  rule over a reference tetrahedron, while the 4th column represent the
%  corresponding weights.
%
% flag: 0: termination with failure, computing rule with possible negative
%          weights and internal nodes
%       1: termination with success, computing rule with positive weights
%          and internal nodes.
%
% iteration_stats: the method perform several iteration and in each one a
%      pointset is analysed so to find the wanted cubature rule that
%      hopefully matches the shifted-tensorial chebyshev moments of total
%      degree "ade"; this variable is a L x 4 matrix in which each rows
%      describes an iteration and is of the form
%
%             [iteration full_set_card inner_set_card residual]
%
%     where:
%
%     * full_set_card: cardinality of the pointset generated in an
%                      hyperrectangle containing the domain;
%
%     * inner_set_card: cardinality of the points of the previous pointset
%                      that are in the polyhedron;
%
%     * residual: how the generated cubature provides moments that are
%                 close to those of the shifted-tensorial chebyshev basis,
%                 based on the infinity norm of their difference.
%
%--------------------------------------------------------------------------
% Routines called by "cubature_polyhedron_tenscheb_02".
%
% 1. rule_search_chebyshev_sub (and its subroutines, main routine is
%    attached below).
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------

% ... troubleshooting ...

if nargin < 5, extraction_type=0; end
if isempty(extraction_type) == 1, extraction_type=0; end
if nargin < 6, set_type=2; end
if isempty(set_type) == 1, set_type=2; end
if nargin < 7, isconvex=0; end
if isempty(isconvex) == 1, isconvex=0; end
if nargin < 8, A=[]; end
if nargin < 9, b=[]; end
if nargin < 10, volume_poly=NaN; end
if size(volume_poly,1) == 0, volume_poly=NaN; end

% ... rule searching ...

[XYZW,flag,iteration_stats,Q,moments1,XYZ_in,V]=rule_search_chebyshev_sub(...
    n,vertices,facets,moments,extraction_type,set_type,isconvex,A,b,...
    volume_poly);

% ... extracting rule by QR method ...

if flag == 0 & (extraction_type ~= 3)
    fprintf('\n \t WARNING: rule with some possible negative weights.  ');
    w=Q'\moments1;
    ind1=find(abs(w)>0); w=w(ind1);
    XYZW=[XYZ_in(ind1,:) w];
    res=norm(V(ind1,:)'*w-moments)/norm(moments);
    
    iteration_stats=[iteration_stats; iteration_stats(end,1:end-1) res];
end









function [XYZW,flag,iteration_stats,Q,moments1,XYZ_in,V]=...
    rule_search_chebyshev_sub(n,vertices,facets,moments,extraction_type,...
    set_type,isconvex,A,b,volume_poly)


%--------------------------------------------------------------------------
% Object:
%
% This routine computes nodes and weights of a cubature formula with degree
% of precision "ade" on a polyhedron defined by "vertices" and "facets".
% The cardinality of the set is at most "(ade+1)*(ade+2)*(ade+3)/6" and the
% nodes are inside the domain.
%
%--------------------------------------------------------------------------
% Input:
%
% n: algebraic degree of precision of the rule (total degree)
%
% vertices: The vertices are stored in "vertices" and is an "m x 3" matrix.
%  The k-th vertex has coordinates "vertices(k,:)".
%
% facets: The variable "facets" can be:
%
% 1. A "cell array":
%    In this case the i-th cell is a vector of positive integers,
%    representing the indices of vertices of the i-th facet, taking into
%    account the variable "vertices". The orientation of the facet is
%    counterclockwise (w.r.t to outward normal).
%
% 2. A "matrix" with dimension "m x 3":
%    In this case the i-th row is a vector of positive integers,
%    representing the indices of vertices of the triangular facet, taking
%    into account the variable "vertices". The orientation of the
%    triangular facet is counterclockwise (w.r.t to outward normal).
%
% extraction_type: algorithm that implements compression:
%    0: Dessole, Marcuzzi, Vianello algorithm (default)
%    1. lsqnonneg
%    2. Slawski's Lawson-Hanson algorithm
%    3. QR compression
%
% set_type: the nodes are taken from certain pointsets:
%    1. tensorial type
%    2. Halton type (default)
%    3. random type
%
% isconvex: if it is known that the polyhedron is convex, then isconvex is
%    equal to 1 otherwise its value is 0 (variable "isconvex" is optional).
%
% A,b: if the polyhedron D is convex then there is a matrix "A" and a
%   vector "b" such that if X=[X(1); X(2); X(3)] belongs to the polyhedron
%   then "A*X<=b" (variables "isconvex" is optional). If they are unknown,
%   set both of them as "[]" (i.e. empty matrix). Both these variables are
%   optional.
%
% volume_poly: approximation of the volume of the domain (it may be unknown
%    and in which case do not set this variable or make it equal to NaN).
%    The variable "volume_poly" is optional.
%
%--------------------------------------------------------------------------
% Output:
%
% XYZW: The variable "XYZW" is a matrix with "m" rows (that is the
%  cardinality of the cubature pointset) and 4 columns.
%  The first 3 columns are the coordinates of the nodes of a
%  rule over a reference tetrahedron, while the 4th column represent the
%  corresponding weights.
%
% flag: 0: termination with failure, computing rule with possible negative
%          weights and internal nodes
%       1: termination with success, computing rule with positive weights
%          and internal nodes.
%
% iteration_stats: the method perform several iteration and in each one a
%      pointset is analysed so to find the wanted cubature rule that
%      hopefully matches the shifted-tensorial chebyshev moments of total
%      degree "ade"; this variable is a L x 4 matrix in which each rows
%      describes an iteration and is of the form
%
%             [iteration full_set_card inner_set_card residual iters]
%
%     where:
%
%     * full_set_card: cardinality of the pointset generated in an
%                      hyperrectangle containing the domain;
%
%     * inner_set_card: cardinality of the points of the previous pointset
%                      that are in the polyhedron;
%
%     * residual: how the generated cubature provides moments that are
%                 close to those of the shifted-tensorial chebyshev basis,
%                 based on the infinity norm of their difference.
%
%     * iters: in the case of Dessole-Marcuzzi-Vianello algorithm or Matlab
%     built-in lsqnonneg it gives the number of inner iterations of the
%     method.
%
%--------------------------------------------------------------------------
% Routines called by "cubature_polyhedron_tenscheb_02".
%
% 1. determine_facets_planes (external)
% 2. test_convexity (attached below)
% 3. chebyshev_moments_polyhedron (and its subroutines, external)
% 4. monomials_moments_convex_polyhedron (and its subroutines, external)
% 5. rule_search_chebyshev_01 (and its subroutines, main routine is
%    attached below).
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------




% ...................... troubleshooting .................................
if nargin < 5, extraction_type=0; end
% 1: meshgrid 2: halton 3. random
if nargin < 6, set_type=2; end
if nargin < 7, isconvex=0; end
if nargin < 8, A=[]; end
if nargin < 9, b=[]; end
if nargin < 10, volume_poly=NaN; end

if isempty(extraction_type), extraction_type=0; end
if isempty(set_type), set_type=2; end
if isempty(isconvex), isconvex=0; A=[]; b=[]; end
if isempty(A), b=[]; end
if isempty(b), A=[]; end

if isconvex == 1
    if isempty(A) | isempty(b)
        [A,b]=determine_facets_planes(vertices,facets);
    end
end

% ......................... initialize ....................................


% Note: on the setting of restol: 10^(-14) seems fine, 10^(-15) is often
% not achieved.
restol=1*10^(-14);

% Active vertices in computations. In case one uses alphashape, some
% vertices may not be used.
ii=unique(facets(:));

bbox(1)=min(vertices(ii,1)); bbox(2)=max(vertices(ii,1));
bbox(3)=min(vertices(ii,2)); bbox(4)=max(vertices(ii,2));
bbox(5)=min(vertices(ii,3)); bbox(6)=max(vertices(ii,3));

XYZ_all=[];   % all points
XYZ_in=[]; % points in the domain
V=[];   % Vandermonde matrix

dim=(n+1)*(n+2)*(n+3)/6; % dimension of polynomial space

% ..... determing pointsets dimensions .....
if isnan(volume_poly) == 1
    rat=10;
    % number of points at each stage
    points_to_generate=rat*8*dim;
else
    rat=(bbox(6)-bbox(5))*(bbox(4)-bbox(3))*(bbox(2)-bbox(1))/volume_poly;
    % number of points at each stage
    points_to_generate=floor(rat*8*dim);
end

fv.vertices=vertices; fv.faces=facets;
k0=0;
iteration_stats=[];
flag=0;

% ......................... exact moments .................................

Lmax=10;
i=0;
in_card=0; card_max=100*dim;

% ... iterations to find cubature rule ...

while (i <= Lmax) & (in_card <= card_max)
    
    % ....................... define local set ............................
    
    % fprintf(' \n \t *');
    
    switch set_type
        case 1
            k=floor(points_to_generate^(1/3));
            x=linspace(bbox(1),bbox(2),k);
            y=linspace(bbox(3),bbox(4),k);
            z=linspace(bbox(5),bbox(6),k);
            [X,Y,Z]=meshgrid(x,y,z);
        case 2
            p=haltonset(3);
            k=points_to_generate;
            
            X=bbox(1)+(bbox(2)-bbox(1))*p(k0+1:k0+k,1);
            Y=bbox(3)+(bbox(4)-bbox(3))*p(k0+1:k0+k,2);
            Z=bbox(5)+(bbox(6)-bbox(5))*p(k0+1:k0+k,3);
            k0=k0+k;
        case 3
            k=points_to_generate;
            X=bbox(1)+(bbox(2)-bbox(1))*rand(k,1);
            Y=bbox(3)+(bbox(4)-bbox(3))*rand(k,1);
            Z=bbox(5)+(bbox(6)-bbox(5))*rand(k,1);
            k0=k0+k;
    end
    
    XYZ_L=[X(:) Y(:) Z(:)];
    
    % .......................... indomain .................................
    
    if isconvex == 0
        inside = inpolyhedron(fv, XYZ_L);
    else
        inside=inconvexpolyhedron(A,b,XYZ_L);
    end
    ind=(inside==1);
    
    
    % .......................... storage ..................................
    
    XYZ_L_in=XYZ_L(ind,:);    % points inside the domain
    
    
    % .......................... find rule ................................
    
    iters=NaN;
    if size(XYZ_L_in,1) < 2*dim
        res=realmax;
    else
        XYZ_in=[XYZ_in; XYZ_L_in];   % all points inside the domain
        VL=chebyshev_vandermonde_matrix(XYZ_L_in,n,bbox);
        V=[V; VL];      % Vandermonde of all points inside the domain
        XYZ_all=[XYZ_all; XYZ_L];      % all points: inside/outside domain
        [Q,R]=qr(V,0);
        
        isn=isnan(R); isnf=max(max(isn));
        
        if norm(imag(R),inf) < 10^(-14)
            R=real(R);
        else
            fprintf('\n \t Warning: complex components in a compression');
            fprintf(' matrix; imag inf norm :%1.1e',norm(imag(R),inf));
            R=real(R);
        end
        
        if norm(imag(Q),inf) < 10^(-14)
            Q=real(Q);
        else
            fprintf('\n \t Warning: complex components in a compression');
            fprintf(' matrix; imag inf norm :%1.1e',norm(imag(Q),inf));
            Q=real(Q);
        end
        
        if norm(imag(moments),inf) < 10^(-14)
            moments=real(moments);
        else
            fprintf('\n \t Warning: complex components in compression');
            fprintf(' moments; imag inf norm :%1.1e',norm(imag(moments),...
                inf));
            moments=real(moments);
        end
        
        moments1=R'\moments;
        
        % .......................... compression ..........................
        
        switch extraction_type
            case 0
                % fprintf('\n \t LHDM');
                [w,resnorm,exitflag, iters] = LHDM(Q', moments1);
                ind1=find(abs(w)>0); w=w(ind1);
            case 1
                % ... w=lsqnonneg(V',moments);
                % fprintf('\n \t lsqnonneg');
                [w,RESNORM,RESIDUAL,EXITFLAG,OUTPUT] =lsqnonneg(Q',...
                    moments1);
                ind1=find(abs(w)>0); w=w(ind1);
                iters=OUTPUT.iterations;
            case 2
                % ... w = lawsonhanson(V', moments);
                %fprintf('\n \t lawsonhanson');
                w = lawsonhanson(Q', moments1);
                ind1=find(abs(w)>0); w=w(ind1);
            otherwise
                % ... w=V'\moments;
                w=Q'\moments1;
                ind1=find(abs(w)>0); w=w(ind1);
        end
        
        
        % ......................... error analysis ........................
        
        res=norm(V(ind1,:)'*w-moments)/norm(moments);
        
        if res>restol | length(w) == 0
            
        else
            XYZW=[XYZ_in(ind1,:) w];
            flag=1;
            iteration_statsL=[i size(XYZ_all,1) size(XYZ_in,1) res iters];
            iteration_stats=[iteration_stats; iteration_statsL];
            return;
        end
        
    end
    
    iteration_statsL=[i size(XYZ_all,1) size(XYZ_in,1) res iters];
    
    iteration_stats=[iteration_stats; iteration_statsL];
    
    
    total_card=size(XYZ_all,1);
    in_card=size(XYZ_in,1);
    
    points_to_generate=next_set_cardinality(total_card,in_card,...
        (4^(i+1))*dim);
    
    i=i+1;
    
end

fprintf('\n \t Compression not completed. Moment error: %1.3e \n',res);
XYZW=[XYZ_in(ind1,:) w];
flag=0;









function V=chebyshev_vandermonde_matrix(XYZ,ade,bbox)

%--------------------------------------------------------------------------
% Object:
% Evaluation of Vandermonde matrix, w.r.t. tensorial Chebyshev matrix in
% the nodes described by the matrix XYZ, whose k-th row has components
% equal to the coordinates of the k-th point of the evaluation.
%--------------------------------------------------------------------------
% Input:
%
% ade: total degree of the basis.
%
% XYZ: The variable "XYZ" is a matrix with "m" rows (that is the
%  cardinality of the cubature pointset) and 3 columns.
%  The latter are the coordinates of the nodes of a
%  rule over a reference tetrahedron, i.e. XYZ(k,:) is a row vector with
%  the coordinates of the k-th point.
%
% bbox: it represents a vector [xm xM ym yM zm zM] so that the tensorial
%  basis is shifted in the hyperrectangle [xm xM] x [ym yM] x [zm zM].
%
%--------------------------------------------------------------------------
% Output:
%
% V: tensorial-Chebyshev Vandermonde matrix. The dimension is
%            number of points x dimension polynomial space
%    hence the k-th row contains the evaluation of the polynomial basis at
%    the k-th point.
%--------------------------------------------------------------------------

XX=XYZ(:,1); YY=XYZ(:,2); ZZ=XYZ(:,3);

a1=bbox(1); b1=bbox(2); A1=(a1+b1)/2; B1=(b1-a1)/2; X=(XX-A1)/B1;
a2=bbox(3); b2=bbox(4); A2=(a2+b2)/2; B2=(b2-a2)/2; Y=(YY-A2)/B2;
a3=bbox(5); b3=bbox(6); A3=(a3+b3)/2; B3=(b3-a3)/2; Z=(ZZ-A3)/B3;

V=[];
for ii=0:ade
    for jj=0:ade-ii
        for kk=0:(ade-ii-jj)
            % VL is a column vector
            VL=cos(ii*acos(X)).*cos(jj*acos(Y)).*cos(kk*acos(Z));
            V=[V VL];
        end
    end
end









function flag=inconvexpolyhedron(A,b,points)

%--------------------------------------------------------------------------
% Object:
%
% Determine if a set of points verifies A*P <= b.
% This is an inpolyhedron routine for a convex polyhedron whose points
% satisfy A*P <=b.
%
%--------------------------------------------------------------------------
% Input:
%
% A: matrix N x 3 describing the facets hyperplanes so that if the domain
% is convex then it satisfies A*[X Y Z]' <= b for each of its point [X Y Z]
%
% b: vector N x 1, determining the domain (if convex)
%
% points: matrix M x 3 (polyhedron vertices coordinates as row vector)
%
%--------------------------------------------------------------------------
% Output:
%
% flag: 0: not convex 1: convex
%
%--------------------------------------------------------------------------

c=A*points';
tol=10^(-14);
points_test= (c <= b+tol);
L=size(A,1);
sum_points_test=sum(points_test,1)-L;
flag=(sum_points_test == 0)';









function next_card=next_set_cardinality(total_card,in_card,card_wanted)

%--------------------------------------------------------------------------
% Object:
%
% Determine the cardinality of a certain pointset, to add at the k-th
% iteration of the rule-search algorithm.
%--------------------------------------------------------------------------
% Input:
%
% total_card: cardinality of all the points generated in a certain
%     hyperrectangle (containing the polyhedron) during the process.
%
% in_card: cardinality of all the points generated in a certain
%     hyperrectangle during the process that are also in the polyhedron.
%
% card_wanted: desired cardinality of nodes that in the next iteration will
%    be in the set.
%
%--------------------------------------------------------------------------
% Output:
%
% next_card: suggested cardinality of nodes to be generated in the
%     hyperrectangle (containing the polyhedron) so that approximatively a
%     number of points equal to "next_card" will be in the polyhedron.
%
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------

% ratio between generated points and points inside the domain

if in_card > 0
    rat=total_card/in_card;
else
    rat=100; % default (no ratio is known, thus many points are suggested)
end

if total_card < card_wanted
    nodes_to_add=card_wanted-total_card;
else
    nodes_to_add=card_wanted; % default
end

next_card=ceil(rat*nodes_to_add);







function [A,b]=determine_facets_planes(vertices,facets)

%--------------------------------------------------------------------------
% Object:
% Determining convex polyhedron equations A*u <= b.
% Each facet suitably determines one of these equations.
%
% A convex polyhedron has points "u=[x; y; z]" that satisfy the equation
% "A*u <= b" for a matrix "A" and a vector "b". We determine such "A" and
% "b", knowing that the k-th facet lies on a plane with equation
% "A(k,:)*u=b(k)".
%
% As procedure we
% 1. compute the equation of each facet;
% 2. determine an other "test" facet;
% 3. see if any the vertices of the test facet satisfy "A(k,:)*u < b(k)"
%    and if not modify "A(k,:)" into "-A(k,:)", "b(k)" into "b(k)".
%--------------------------------------------------------------------------
% Input:
%
% vertices: The vertices are stored in "vertices" and is an "m x 3" matrix.
%  The k-th vertex has coordinates "vertices(k,:)".
%
% facets: The variable "facets" can be:
%
% 1. A "cell array":
%    In this case the i-th cell is a vector of positive integers,
%    representing the indices of vertices of the i-th facet, taking into
%    account the variable "vertices". The orientation of the facet is
%    counterclockwise (w.r.t to outward normal).
%
% 2. A "matrix" with dimension "m x 3":
%    In this case the i-th row is a vector of positive integers,
%    representing the indices of vertices of the triangular facet, taking
%    into account the variable "vertices". The orientation of the
%    triangular facet is counterclockwise (w.r.t to outward normal).
%--------------------------------------------------------------------------
% Output:
%
% A,b: matrix and vector so that the k-th facet is contained in the plane
%     whose points X=(x(1); x(2); x(3)) satisfy the equation
%
%                              A(k,:)*X <= b(k)
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------

tf = iscell(facets);
number_of_facets=length(facets);
A=[]; b=[];

for k=1:number_of_facets
    
    % Row vector of indices of vertices describing the k-th facet.
    kt=setdiff(1:number_of_facets,k);
    
    if tf == 1 % facets are stored in cells
        facetL=facets{k}; facet_test_L=facets{kt};
    else % facets are stored in m x 3 matrix (triangular facets)
        facetL=facets(k,:);
        facet_test_L=facets(kt,:);
    end
    
    % Vertices of the facet ordered counterclockwise w.r.t. the outward
    % normal.
    XV=vertices(facetL',1); YV=vertices(facetL',2); ZV=vertices(facetL',3);
    
    % ..... Computing plane equation .....
    
    P1=[XV(1) YV(1) ZV(1)]; P1=P1';
    P2=[XV(2) YV(2) ZV(2)]; P2=P2';
    P3=[XV(3) YV(3) ZV(3)]; P3=P3';
    
    % Equation of the plane containing the facets.
    % See e.g.:
    % https://math.stackexchange.com/questions/2686606/...
    % equation-of-a-plane-passing-through-3-points
    %
    % Here "a" is a column vector and "b" is a scalar.
    
    aL=cross(P2-P1,P3-P1); bL=P1'*aL;
    
    XV_test=vertices(facet_test_L',1);
    YV_test=vertices(facet_test_L',2);
    ZV_test=vertices(facet_test_L',3);
    
    % Warning:
    % This check is a little dangerous due to possible machine errors.
    % Must be improved in next releases.
    residual_test_vertices0=bL-aL'*[XV'; YV'; ZV'];
    
    residual_test_vertices=bL-aL'*[XV_test'; YV_test'; ZV_test'];
    [max_distance_test_vertex,index]=max(abs(residual_test_vertices));
    
    convex_check=residual_test_vertices(index);
    
    
    if convex_check < 0
        % warning('Modifying plane equation.');
        aL=-aL; bL=-bL;
    else
        if convex_check == 0
            warning('Troubles in determining polyhedron equation Ax<=b.');
        end
    end
    
    A=[A; aL']; b=[b; bL];
    
end










function moments=chebyshev_moments_polyhedron(vertices,facets,ade,...
    normals)

%--------------------------------------------------------------------------
% Object:
%
% Compute the moments "moments" of the chebyshev basis
%                         T_i(x)*T_j(y)*T_k(z)
% with i+j+k <= ade, over the polyhedron with vertices "vertices", whose
% polygonal facets are described by "facets".
%--------------------------------------------------------------------------
% Input: 
% vertices: The vertices are stored in "vertices" and is an "m x 3" matrix.
% The k-th vertex has coordinates "vertices(k,:)".
%
% facets: The variable "facets" can be:

% 1. A "cell array":
%    In this case the i-th cell is a vector of positive integers,
%    representing the indices of vertices of the i-th facet, taking into
%    account the variable "vertices". The orientation of the facet is
%    counterclockwise (w.r.t to outward normal).
%
% 2. A "matrix" with dimension "m x 3":
%    In this case the i-th row is a vector of positive integers,
%    representing the indices of vertices of the triangular facet, taking
%    into account the variable "vertices". The orientation of the
%    triangular facet is counterclockwise (w.r.t to outward normal).
%
% If facets normals are available they are stored in the matrix "normals".
% The k-th row contains the normal of the k-th facet.
%
% ade: algebraic (total) degree of precision of the rule.
%
% normals: normals to the facets.
%--------------------------------------------------------------------------
% Output:
% moments: moments of the chebyshev basis
%                         T_i(x)*T_j(y)*T_k(z)
% with i+j+k <= ade, over the polyhedron with vertices "vertices", whose
% polygonal facets are described by "facets".
%--------------------------------------------------------------------------
% Reference:
% [1] E.B. Chin, J.B. Lasserre, N. Sukumar:
% "Numerical integration of homogeneous function on convex and nonconvex
% polygons and polyhedra".
% Computational Mechanics, Vol. 56, No. 6, pp 967?981.
% https://doi.org/10.1007/s00466-015-1213-7
%
% See also references therein.
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
% Chebyshev indices.
%
% Since we have to compute the cubature of a basis of tensorial Chebyshev
% polynomials of total degree up to "ade", we have to provide an order to
% the Chebyshev indices.
%
% More in detail, they are of the form
%                  "(T_i(x))*(T_j(y))*(T_k(z))"
% and we order such triples (i,j,k).
%--------------------------------------------------------------------------

chebyshev_indices=[];
for ii=0:ade
    for jj=0:ade-ii
        for kk=0:(ade-ii-jj)
            chebyshev_indices=[chebyshev_indices; ii jj kk];
        end
    end
end

tf = iscell(facets);
number_of_facets=length(facets);



%--------------------------------------------------------------------------
% Computing hyperrectangle with sides parallel to the x, y, z axis.
%
% Active vertices in computations. In case one uses alphashape, some
% vertices may not be used.
%--------------------------------------------------------------------------

ii=unique(facets(:));

bbox(1)=min(vertices(ii,1)); bbox(2)=max(vertices(ii,1));
bbox(3)=min(vertices(ii,2)); bbox(4)=max(vertices(ii,2));
bbox(5)=min(vertices(ii,3)); bbox(6)=max(vertices(ii,3));



%--------------------------------------------------------------------------
% Computing moments "chebyshev_moms" on polyhedron facets.
% The moments are ordered following "chebyshev_indices" to determine the
% m-th polynomial to be considered.
% The matrix "chebyshev_moms" has "number_of_facets" columns and the
% dimension of the polynomial space of total degree "ade" as number of rows.
%
% In other words "chebyshev_moms(ii,jj)" is the moment of the "ii"-th
% chebyshev basis on the "jj"-th facet.
%--------------------------------------------------------------------------


chebyshev_moms=[];
for k=1:number_of_facets
    
    % Row vector of indices of vertices describing the k-th facet
    if tf == 1 % facets are stored in cells
        facetL=facets{k};
    else % facets are stored in m x 3 matrix (triangular facets)
        facetL=facets(k,:);
    end
    
    % Vertices of the facet ordered counterclockwise w.r.t. the outward
    % normal.
    XV=vertices(facetL',1);
    YV=vertices(facetL',2);
    ZV=vertices(facetL',3);
    
    % Cubature rule over facet. Storing in "n" the normal to the facet.
    if nargin < 4
        [XYZW,n]=cubature_facet_formula(XV,YV,ZV,ade+1);
    else
        if size(normals,1) == 0
            [XYZW,n]=cubature_facet_formula(XV,YV,ZV,ade+1);
        else
            [XYZW]=cubature_facet_formula(XV,YV,ZV,ade+1);
            n=normals(k,:);
        end
    end
    
    nodes=XYZW(:,1:3); weights=XYZW(:,4);
    chebyshev_moms_facet=cubature_tens_chebyshev_facet(nodes,weights,...
        chebyshev_indices, bbox);
    
    
    
    % Note:
    % The factor n(1) comes out from the divergence theorem applied to the
    % vector field V=(f,0,0) where f is the x-primitive of the
    % tensorial-Chebyshev polynomial of which one wants to compute the
    % moment.
    
    
    chebyshev_moms=[chebyshev_moms n(1)*chebyshev_moms_facet];
    
end




%--------------------------------------------------------------------------
% Computing moments "moments" on polyhedron by summing contributions on the
% facets.
% See reference [1], formula (7).
%--------------------------------------------------------------------------

moments=sum(chebyshev_moms,2);








function chebyshev_moms=cubature_tens_chebyshev_facet(nodes,weights,...
    chebyshev_indices,bbox)

% Compute the cubature on the facet, having cubature nodes/weights of an
% algebraic rule with nodes "nodes" and weights "weights" of degree of
% precision "ade", of all the chebyshevs having total degree inferior or
% equal to "ade".
% As output we provide these moments in the column vector "chebyshev_moms",
% ordered as in "chebyshev_indices", with total degree as
% "chebyshev_total_degrees".

X=nodes(:,1); Y=nodes(:,2); Z=nodes(:,3);

% Map to reference cube [-1,1] x [-1,1] x [-1,1]
a1=bbox(1); b1=bbox(2); A1=(a1+b1)/2; B1=(b1-a1)/2; XN=(X-A1)/B1;
a2=bbox(3); b2=bbox(4); A2=(a2+b2)/2; B2=(b2-a2)/2; YN=(Y-A2)/B2;
a3=bbox(5); b3=bbox(6); A3=(a3+b3)/2; B3=(b3-a3)/2; ZN=(Z-A3)/B3;

chebyshev_moms=[];

for index=1:size(chebyshev_indices,1)
    
    i=chebyshev_indices(index,1);
    j=chebyshev_indices(index,2);
    k=chebyshev_indices(index,3);
    
    % Tensorial Chebyshev polynomials x-primitive.
    switch i
        case 0
            f=@(x,y,z) B1*x.*cos(j*acos(y)).*cos(k*acos(z));
        case 1
            f=@(x,y,z) B1*((x.^2)/2).*cos(j*acos(y)).*cos(k*acos(z));
        otherwise
            f=@(x,y,z) B1*((i*cos((i+1)*acos(x))/(i^2-1)- ...
                x.*cos((i)*acos(x))/(i-1)).*cos(j*acos(y)).* ...
                cos(k*acos(z)));
    end
    
    FXYZN=feval(f,XN,YN,ZN);
    chebyshev_momsL=weights'*FXYZN;
    chebyshev_moms=[chebyshev_moms; chebyshev_momsL];
    
end









function [XYZW,n]=cubature_facet_formula(XV,YV,ZV,ade)

%--------------------------------------------------------------------------
% OBJECT:
% This routine computes a cubature rule of algebraic degree of exactness
% equal to "ade" on the polygon in 3D whose vertices are (XV,YV,ZV).
% As output we have a matrix "XYZW" with "m" rows and "4" columns.
% The nodes are the rows of "XYZW(:,1:3)", while the weights are
% "XYZW(:,4)".
%
% n is the outer normal vector to the facet of the polyhedron, supposed
% that the vertices are oriented counterclockwise in its respect.
%
% Note: differently from classical codes on polygons, the first and
% last vertex must NOT be equal, i.e. the points [XV(1) YV(1) ZV(1)] and
% [XV(end) YV(end) ZV(end)] must not be equal.
%--------------------------------------------------------------------------



% Rototranslation map to a certain 2D polygon
[xv0,yv0,zv0,R10,RFM]=maptopolygon2(XV,YV,ZV);

% Cubature over a certain 2D polygon.
[xyw,xvc,yvc,pgon,tri]=polygauss_2018(ade,xv0,yv0);

xv_cub=xyw(:,1); yv_cub=xyw(:,2); zw_cub=zv0(1)*ones(size(yv_cub));

% Map it back to obtain a cubature rule over the 3D polygon.
[XV_CUB,YV_CUB,ZV_CUB,n]=maptopolygon3(xv_cub,yv_cub,zw_cub,R10,RFM);
WV_CUB=xyw(:,3);

% Assembly the rule as a certain matrix.
XYZW=[XV_CUB YV_CUB ZV_CUB WV_CUB];









function area_triangle=area_triangular_facet(XV,YV,ZV)

% AREA OF A TRIANGLE IN 3D VIA ERONE'S RULE.

A=[XV(1) YV(1) ZV(1)];
B=[XV(2) YV(2) ZV(2)];
C=[XV(3) YV(3) ZV(3)];

c=norm(B-A); a=norm(C-B); b=norm(C-A);
p=(a+b+c)/2;

% Erone formula.
area_triangle=sqrt(p*(p-a)*(p-b)*(p-c));









function [xv,yv,zv,R10,RFM,abcd]=maptopolygon2(XV,YV,ZV)

%--------------------------------------------------------------------------
% Object:
% This routine takes a polygon in the 3D space whose vertices are (XV,YV,ZV)
% and maps it to a flat polygon with vertices (xv,yv,0) by a "rigid
% motion".
% The vector "zv" and the rotation matrix "R10" allow to pull back the
% points.
%--------------------------------------------------------------------------
% Input:
%
% XV,YV,ZV: coordinates of the vertices of a polygon in the three
% dimensional space (e.g. the vertices of a facet of a tetrahedra); last
% vertex may not be repeated.
%--------------------------------------------------------------------------
% Output:
%
% xv,yv: coordinates of the vertices of a polygon in the two
% dimensional space that is obtained by a "rigid motion" of a polygon with
% vertices (XV,YV,ZV); last vertex may not be repeated.
%
% zv, R10: the polygon P with coordinates (xv,yv,0) is
%          (R10*[XV'; YV'; ZV'])'-[0V 0V ZV]
%     where 0V is the null vector with the same dimension of ZV.
%
% R10,RFM: rotation and reflection matrices involved in the process.
%     It is "P=R10'*(RFM*P2+PS)" where "PS" is a certain shift, "P2" a
%     the plane xy and "P" is a point in the 3D plane.
%
% abcd: vector whose components [a b c d] determine the equation of the
%     plane "a*x+b*y+c*z+d=0" containing the vertices of the polygon.
%
%--------------------------------------------------------------------------

% Determine the needed rigid motion. We map A in 0, B in a point with
% coordinates (dist(A,B),0,0) and the versor of the normal into (0,0,1).

A=[XV(1) YV(1) ZV(1)]; A=A';
B=[XV(2) YV(2) ZV(2)]; B=B';
C=[XV(3) YV(3) ZV(3)]; C=C';

% Equation of the plane containing the facets.
% See e.g.:
% https://math.stackexchange.com/questions/2686606/equation-of-a-plane-passing-through-3-points
n=cross(B-A,C-A);
a=n(1); b=n(2); c=n(3);
d=-A'*n;

abcd=[a b c d];

% Rotation to the xy plane.
% https://math.stackexchange.com/questions/1167717/transform-a-plane-to-the-xy-plane
% Consider correction due to user "paterry".

if (a^2+b^2) > 0 % rotation is required
    den=sqrt(a^2+b^2+c^2);
    cos_theta=c/den;
    sin_theta=sqrt(a^2+b^2)/den;
    
    den2=sqrt(a^2+b^2);
    u1=b/den2;
    u2=-a/den2;
    
    term_cos=(1-cos_theta);
    
    R_row1=[cos_theta+u1^2*term_cos  u1*u2*term_cos u2*sin_theta];
    R_row2=[u1*u2*term_cos cos_theta+u2^2*term_cos  -u1*sin_theta];
    R_row3=[-u2*sin_theta u1*sin_theta cos_theta];
    R10=[R_row1; R_row2; R_row3];
else % no rotation is required
    R10=eye(3);
end

xyv=R10*[XV'; YV'; ZV']; xyv=xyv';

xv=xyv(:,1);
yv=xyv(:,2);
zv=xyv(:,3);

% Check orientation polygon
[orient,signed_area] = polyorient(xv,yv);

if orient == 0
    % if orientation is clockwise, we reflect the 2D polygon w.r.t. the
    % x-axis so to map it into an equivalent one with the right orientation
    RFM=diag([1 -1 1]);
    yv=-yv;
    [orient,signed_area] = polyorient(xv,yv);
else
    RFM=eye(3);
end

if orient == 0
    warning('Some issues with orientation');
end






 






function [orient,signed_area] = polyorient(x,y)
%--------------------------------------------------------------------------
%   POLYORIENT Orientation of polygon
%   Returns the orientation and signed area of a 2D polygon
%
%   Syntax:
%      [orient,signed_area] = POLYORIENT(X,Y)
%
%   Inputs:
%      X, Y   Vectors with the polygon vertices
%
%   Outputs:
%      orient   Polygon orientation. 1 if the orientation is 
%            counter-clockwise (direct), 0 otherwise
%      signed_area    Signed area of the polygon, negative if orientation 
%               is not direct
%--------------------------------------------------------------------------
%   Examples:
%      x1 = [0 0 1 1]; y1 = [1 2 2 1];
%      x2 = [0 0 1 1]; y2 = [1 0 0 1];
%      x3 = [x1 x2];   y3 = [y1 y2];
%
%      [o1,a1] = polyorient(x1,y1) % 0, -1
%      [o2,a2] = polyorient(x2,y2) % 1,  1
%      [o3,a3] = polyorient(x3,y3) % 0,  0
%--------------------------------------------------------------------------
% Author: 
% MMA 21-4-2006, mma@odyle.net
%--------------------------------------------------------------------------

signed_area=0.5* sum(x.*y([2:end 1]) - y.*x([2:end 1]));
orient = signed_area > 0;









function [XV,YV,ZV,n]=maptopolygon3(xv,yv,zv,R10,RFM)

%--------------------------------------------------------------------------
% Object:
% This routine takes a planar polygon in the 3D space whose vertices are 
% (XV,YV,ZV) and maps it to a flat polygon with vertices (xv,yv,0) by a 
% "rigid motion".
% The vector "zv" and the rotation matrix "R10" allow to map back the
% points.
%--------------------------------------------------------------------------
% Input:
%
% xv, yv: coordinates of the vertices of a polygon in the two dimensional
%   space that is obtained by a "rigid motion" of a polygon with
%   vertices (XV,YV,ZV); last vertex may not be repeated; they are column 
%   vectors.
%
% R10,RFM: rotation and reflection matrices involved in the process.
%     It is "P=R10'*(RFM*P2+PS)" where "PS" is a certain shift, "P2" a
%     the plane xy and "P" is a point in the 3D plane.
%
%--------------------------------------------------------------------------
% Output:
%
% XV,YV,ZV: coordinates of the vertices of a polygon in the three
% dimensional space (e.g. the vertices of a facet of a tetrahedra); last
% vertex may not be repeated; they are column vectors.
%
% n: versor of outer normal to the plane containing [XV YV ZV].
%
%--------------------------------------------------------------------------
% Data:
% Revised: February 16, 2021.
%--------------------------------------------------------------------------

% k-th column determines the k-th point
P2=[xv'; yv'; zeros(size(zv'))];
PS=[zeros(size(zv')); zeros(size(zv')); zv' ];

P=R10'*(RFM*P2+PS);

XV=(P(1,:))';
YV=(P(2,:))';
ZV=(P(3,:))';

if nargout >= 4
    ee=[0;0;1];
    n=R10'*(RFM*ee); % outer normal to the facet
    n=n/norm(n); % versor of outer normal
end











