
function [moments,A,b]=monomials_moments_convex_polyhedron(vertices,...
    facets,ade,A,b)

%--------------------------------------------------------------------------
% Object:
%
% Compute the moments "monomial_moments" of the monomial basis
%                         x^i * y^j * z^k
% 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).
%
% ade: computation of all the moments of monomial basis (total degree at
%    most "ade").
%
% A,b: each point X=[x(1); x(2); x(3)] satifies A*X<=b (this input is not
%    mandatory).
% 
%--------------------------------------------------------------------------
% Output:
%
% moments: moments of the monomial basis, in lexicographic order, e.g. for
%   "ade=2" the monomial exponents "i", "j", "k" are
%
%      0     0     0
%      0     0     1
%      0     0     2
%      0     1     0
%      0     1     1
%      0     2     0
%      1     0     0
%      1     0     1
%      1     1     0
%      2     0     0
%
% A,b: respectively matrix and vector determining the convex polyhedron D,
%   i.e. X=[x; y; z] belongs to the domain D if and only if "A*X<=b".
%
%--------------------------------------------------------------------------
% Subroutines called by "monomials_moments_convex_polyhedron":
%
% 1. cubature_monomials_facet (attached)
% 2. cubature_facet_formula (attached)
% 3. cubature_facet_formula (attached,somehow calls:
%    * area_triangular_facet (attached)
%    * cubature_triangle
%    * maptopolygon2 (and its subroutines)
%    * maptopolygon3
%    * polygauss_2018
% 4. determine_facets_planes 
% 
%--------------------------------------------------------------------------
% 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.
%--------------------------------------------------------------------------


if nargin < 4
    A=[]; n=[];
end




%--------------------------------------------------------------------------
% Monomial exponentials.
%
% Since we have to compute the cubature of a basis of monomials of total
% degree up to "ade", we have to provide an order to the monomials
% exponents.
%
% More in detail, they are of the form "(x^i)*(y^j)*(z^k)" and we order
% such triples (i,j,k).
%--------------------------------------------------------------------------

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

% Determining convex polyhedron equations A*u <= b.
if size(A,1) == 0
[A,b]=determine_facets_planes(vertices,facets);
end


%--------------------------------------------------------------------------
% Computing moments "monomial_moms" on polyhedron facets.
% The moments are ordered following "monomial_exponents" to determine the
% m-th monomial to be considered.
% The matrix "monomial_moms" has "number_of_facets" columns and the
% dimension of the polynomial space of total degree "ade" as number of rows.
%
% In other words "monomial_moms(ii,jj)" is the moment of the "ii"-th
% monomial basis on the "jj"-th facet.
%
% Technical note: the routine below can be made faster to compute the 
% overall moments on the polyhedron by not computing moments when on the 
% facet we have "b=0", in view of formula (7) in
% E.B. Chin, J.B. Lasserre, N. Sukumar reference.
%--------------------------------------------------------------------------

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

monomial_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.
    XYZW=cubature_facet_formula(XV,YV,ZV,ade);
    
    nodes=XYZW(:,1:3); weights=XYZW(:,4);
    monomial_moms_facet=cubature_monomials_facet(nodes,weights,ade,...
        monomial_exponents);
    monomial_moms=[monomial_moms monomial_moms_facet];
end




%--------------------------------------------------------------------------
% Computing moments "monomial_moms" on polyhedron.
% See reference [1], formula (7).
%--------------------------------------------------------------------------

p=3; q=sum(monomial_exponents,2); 
outer_coeff=1./(p+q);

den_ic=sqrt(A(:,1).^2+A(:,2).^2+A(:,3).^2);
inner_coeff=b./den_ic;

monomial_moms(1,:);
moments=outer_coeff.*(monomial_moms*inner_coeff);









function monomial_moms=cubature_monomials_facet(nodes,weights,ade,...
    monomial_exponents)

% 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 monomials having total degree inferior or
% equal to "ade".
% As output we provide these moments in the column vector "monomial_moms",
% ordered as in "monomial_exponents", with total degree as
% "monomial_total_degrees".

XN=nodes(:,1); YN=nodes(:,2); ZN=nodes(:,3);
monomial_moms=[];
for index=1:size(monomial_exponents,1)
    
    ii=monomial_exponents(index,1);
    jj=monomial_exponents(index,2);
    kk=monomial_exponents(index,3);
    
    FXYZN=(XN.^ii).*(YN.^jj).*(ZN.^kk);
    monomial_momsL=weights'*FXYZN;
    monomial_moms=[monomial_moms; monomial_momsL];
    
end







function XYZW=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)".
%
% 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.
%--------------------------------------------------------------------------

if length(XV) == 3
    
    % THE POLYGON IS A TRIANGLE.
    
    [xyw,xyw_bar]=cubature_triangle(ade);
    L=xyw_bar(:,1:3);
    XYZ_CUB=L*[XV YV ZV];
    area_triangle=area_triangular_facet(XV,YV,ZV);
    WV_CUB=area_triangle*xyw_bar(:,4);
    XYZW=[XYZ_CUB WV_CUB];
    
else
    
    % THE POLYGON IS NOT A TRIANGLE.
    
    % 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]=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];
    
end







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));


