function dec = appendUnifiedPartition(dec, seIndexVector, zerotol, ...
    isValid, isConfined)
% add unified partition, calculate simplex partitions and perform book
%   keeping. Up to here, I know the extreme rays of the unified
%   partition.
% The strategy is:
%   (1) Triangulate the vrep of the unified partition
%   (2) Generate the simplex partitions (producing hreps from vreps is
%       simple)
%   (3) Filter the outer facets of the unified partition from the
%       simplex partitions.
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% 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 <http://www.gnu.org/licenses/>.

%! Note: Why was the conversion vrep->hrep based on the conversion
%           function?
% Coming from the vrep, that could be given by:
%     struct('V', [dec.seScaledToClose(seIndexVector, :); ...
%                  zeros(1, dec.nC)]);
%   is not good. Testing cs-decomposition throws errors for the
%   conversion afterwards, which is:
%     partition = convertVrepToHrep(vrep, zerotol);
% It was implemented like above, but changed for the given reasons.
%/

%% initialization
if ~isfield(dec, 'unifiedPartition') || ~exist('seIndexVector', 'var')
    dec.simplexPartition = cell(0);
    dec.unifiedPartition = cell(0);
    dec.simplexToUnifiedVector = [];
    dec.simplexIsValidVector = [];
    dec.simplexIsConfinedVector = [];
    dec.unifiedToSimplexVector = [];
    dec.unifiedIsValidVector = [];    
    dec.unifiedIsConfinedVector = [];
    dec.seIncidenceMatrix = eye(dec.nS, dec.nS);
    return
end

thisUnifiedIndex = length(dec.unifiedPartition) + 1;

% triangulation of unified partition
seSimplexIndexMatrix = computeTriangulationOfCone(...
    dec.seScaledToClose(seIndexVector, :));

% prepare matrix of facet normals from all simplex partitions
simplexNormalVectorMatrix = [];
simplexVolumeSum = 0;
    
% adding simplex partitions
for iSimplex = 1:size(seSimplexIndexMatrix, 1)
    
    % that's the new index for the simplex partition
    thisSimplexIndex = length(dec.simplexPartition) + 1;
    
    % thats the matrix of structuring elements
    thisSeMatrix = dec.seScaledToClose(seIndexVector(...
        seSimplexIndexMatrix(iSimplex, :)), :);
    
    % sometimes, it happens that the triangulation is not full. We
    %   ignore these components.
    if rank(thisSeMatrix, zerotol) < dec.nC
        continue
    end
    
    % For simplex partitions, the hrep can be generated directly by inversion
    %   of the matrix of structuring elements. For a unified partition, this is
    %   not possible. How does this work?
    % The points of the cone is described by:
    %     {hc = se' * lam | lam >= 0}
    %   where the equation can be converted to:
    %     lam = inv(se') * hc
    %   to obtain the corresponding H-representation:
    %     {hc | -inv(se') * hc <= 0}
    
    % generate hrep of partition
    dec.simplexPartition{thisSimplexIndex}.A = -1*inv(thisSeMatrix');
    dec.simplexPartition{thisSimplexIndex}.h = zeros([dec.nC, 1]);
    dec.simplexPartition{thisSimplexIndex} = normalizeHrep(...
        dec.simplexPartition{thisSimplexIndex});
    
    % compose the matrix of all normal vectors
    simplexNormalVectorMatrix = [simplexNormalVectorMatrix; ...
        dec.simplexPartition{thisSimplexIndex}.A];
    
    % add structuring element indices
    dec.simplexPartition{thisSimplexIndex}.seIndexVector = ...
        sort(seIndexVector(seSimplexIndexMatrix(iSimplex, :)));
    
    % obtain scaled volume
    dec.simplexPartition{thisSimplexIndex}.volumeScaled = ...
        computeVolumeVrep(...
        struct('V', [thisSeMatrix; zeros([1, dec.nC])]) ...
        ) / dec.confinementCone.volume;
    simplexVolumeSum = simplexVolumeSum + ...
        dec.simplexPartition{thisSimplexIndex}.volumeScaled;
    
    % and we need to add this index
    dec.simplexToUnifiedVector(thisSimplexIndex) = thisUnifiedIndex;
    dec.simplexIsValidVector(thisSimplexIndex) = isValid;
    dec.simplexIsConfinedVector(thisSimplexIndex) = isConfined;
    
    % and we should add the mapping
    dec.simplexPartition{thisSimplexIndex}.mappingHcToLam = inv(dec.se(...
        dec.simplexPartition{thisSimplexIndex}.seIndexVector, : ...
        )');
end

% The simplex partitions were generated from a triangulation from which
%   it is assumed that this was complete. As the unified partition is
%   the union of all simplex partitions, we only have to filter for
%   facets normals that have no opposing normal vector. These must be
%   the outer normal vectors of the union.

% generate unified hrep from vrep
dec.unifiedPartition{thisUnifiedIndex}.A = eliminateRows(...
    simplexNormalVectorMatrix, -1*simplexNormalVectorMatrix, zerotol);
dec.unifiedPartition{thisUnifiedIndex}.h = zeros(...
    [size(dec.unifiedPartition{thisUnifiedIndex}.A, 1), 1]);

% add structuring element indices
dec.unifiedPartition{thisUnifiedIndex}.seIndexVector = sort(seIndexVector);

% add mapping to simplex partition (simply the last simplex mapping entry)
dec.unifiedToSimplexVector(thisUnifiedIndex) = length(dec.simplexToUnifiedVector);
dec.unifiedIsValidVector(thisUnifiedIndex) = isValid;
dec.unifiedIsConfinedVector(thisUnifiedIndex) = isConfined;

% obtain scaled volume
dec.unifiedPartition{thisUnifiedIndex}.volumeScaled = ...
    computeVolumeVrep(...
    struct('V', [...
    dec.seScaledToClose(seIndexVector, :); ...
    zeros([1, dec.nC])...
    ]) ...
    ) / dec.confinementCone.volume;

% add incidence information
for iSeIndex = 1:length(seIndexVector)
    for jSeIndex = 1:length(seIndexVector)
        dec.seIncidenceMatrix(seIndexVector(iSeIndex), ...
            seIndexVector(jSeIndex)) = 1;
    end
end

end