function dec = addConfinementPartitionData(dec, partitionData, zerotol)
% add confined partitions
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 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/>.

%% construct and accumulate SEs
% each partition contains only the open hrep, we need to construct the
%   vrep, extract the ray elements and append them as SEs
dec.seScaledToClose = zeros(0, dec.nC);
for iPart = 1:length(partitionData)
%     iPart
    % skip unconfined partitions
    if ~partitionData(iPart).isConfined,continue;end;
    
    % get current structuring elements
    thisHrepClosed = closeCone(partitionData(iPart).hrep, dec.centerRay);
    thisVrepClosed = convertHrepToVrep(thisHrepClosed);
    thisVrep = openCone(thisVrepClosed, zerotol);
%     thisSeMatrix = scaleToReal(thisVrep.V, dec, zerotol)
        
    % filter known SEs
    [thisSeMatrix thisSeIndexVector] = eliminateRows(thisVrep.V, ...
        dec.seScaledToClose, zerotol);
    % append new SEs
    thisSeIndexVector = [thisSeIndexVector, size(dec.seScaledToClose, 1) + ...
        (1:size(thisSeMatrix, 1))];
    dec.seScaledToClose = [dec.seScaledToClose; thisSeMatrix];
    % store result
    partitionData(iPart).seIndexVector = sort(thisSeIndexVector);
end
% scale all structuring elements to physical
dec.se = scaleToReal(dec.seScaledToClose, dec, zerotol);

%% recalculate center ray
dec.centerRay = mean(dec.se, 1)';
dec.centerRay = dec.centerRay / norm(dec.centerRay);
% scale SEs
dec.seScaledToClose = scaleToClose(dec.se, dec);
% volume of validity cone is, hence, invalid - instead of recalculating it,
%   it is deleted
dec.validityCone = rmfield(dec.validityCone, 'volume');
% se indices
dec.validityCone.seIndexVector = find(any(abs(dec.validityCone.A * ...
    dec.seScaledToClose') < zerotol, 1));

%% construct confinement cone
% all structuring elements span the confinement cone, even though it is not
%   a proper positive basis.
vrepClosedCone.V = closeCone(dec.seScaledToClose);
hrepClosedCone = convertVrepToHrep(vrepClosedCone, zerotol);
dec.confinementCone = openCone(hrepClosedCone);
% add SE indices
dec.confinementCone.seIndexVector = find(any(abs(dec.confinementCone.A * ...
    dec.seScaledToClose') < zerotol, 1));
% add cone volume
dec.confinementCone.volume = computeVolumeVrep(vrepClosedCone);

%% add unified & simplex partitions
% we will track simultaneaously the total volume of all unified partitions.
%   Which must sum up to the confinement cone volume for non-odd
%   representations
% initialize
dec = appendUnifiedPartition(dec);
unifiedVolumeSum = 0;
% execute
for iPart = 1:length(partitionData)
    % skip if not confined
    if ~partitionData(iPart).isConfined
        continue
    end
    
    % add partition as usual   
    dec = appendUnifiedPartition(dec, partitionData(iPart).seIndexVector, ...
        zerotol, partitionData(iPart).isValid, 1);
    
    dec = appendPartitionDataForEmbedding(dec, partitionData(iPart), zerotol);
    
    % accumulate unified partition volume
    unifiedVolumeSum = unifiedVolumeSum + ...
        dec.unifiedPartition{end}.volumeScaled;
    
end
dec.nPU = length(dec.unifiedPartition);
dec.nPS = length(dec.simplexPartition);
dec.nS = size(dec.se, 1);

%% check oddity
% the confinement cone for odd representations is not convex. Here, the
%   scaled volumes of the unified partitions don't sum up to one.
if abs(unifiedVolumeSum - 1) <= zerotol
    dec.isOrdinary = 1;
    return
else
    dec.isOrdinary = 0;
    % if the result is singular, we must add some more partitions
end

%% add unconfined partitions
for iPart = 1:length(partitionData)
    % skip if confined
    if partitionData(iPart).isConfined, continue;end;
    
    % compute SEs (physically scaled)
    thisHrepClosed = closeCone(struct(...
        'A', [dec.confinementCone.A; partitionData(iPart).hrep.A], ...
        'h', [dec.confinementCone.h; partitionData(iPart).hrep.h]), ...
        dec.centerRay);
    thisVrepClosed = convertHrepToVrep(thisHrepClosed);
    thisVrep = openCone(thisVrepClosed, zerotol);
    thisSeMatrix = scaleToReal(thisVrep.V, dec, zerotol);
    
    % as we intersect a partition with the confinement cone the result can
    %   be a lower dimensional entity at the boundary of the confinement
    %   cone.
    if computeVolumeVrep(thisVrepClosed)/dec.confinementCone.volume < zerotol
        continue
    end
    
    % We expect that a confined partition is either fully contained in the
    %   confinement cone or it is fully outside of the confinement cone.
    %   Hence, we check if the extreme rays reside fully outside or fully
    %   inside:
    distanceMatrix = dec.confinementCone.A * thisSeMatrix';
    if all(distanceMatrix(:) > -zerotol)
        % partition completely outside
        continue
    elseif all(distanceMatrix(:) < zerotol)
        % partition completely inside, continue below
    else
        error('elk:decomposition:numericalError', ['I''m confused since a ' ...
            'partition should either fully reside inside the confinement ' ...
            'cone or fully outside. Most likely, the numerics fail, here. ' ...
            'This should not happen.']);
    end
        
    % filter known SEs
    [thisSeMatrix thisSeIndexVector] = eliminateRows(thisSeMatrix, ...
        dec.se, zerotol);
    % append new SEs
    thisSeIndexVector = [thisSeIndexVector, size(dec.se, 1) + ...
        (1:size(thisSeMatrix, 1))];
    dec.se = [dec.se; thisSeMatrix];
    dec.seScaledToClose = [dec.seScaledToClose; ...
                           scaleToClose(thisSeMatrix, dec)];
    
    % add partition as usual
    dec = appendUnifiedPartition(dec, thisSeIndexVector, zerotol, ...
        partitionData(iPart).isValid, 0);
    dec = appendPartitionDataForEmbedding(dec, partitionData(iPart), zerotol);
      
    % accumulate unified partition volume
    unifiedVolumeSum = unifiedVolumeSum + ...
        dec.unifiedPartition{end}.volumeScaled;
      
end

% recompute some values:
dec.nPU = length(dec.unifiedPartition);
dec.nPS = length(dec.simplexPartition);
dec.nS = size(dec.se, 1);

% if abs(unifiedVolumeSum - 1) > zerotol
%     error('elk:decomposition:numericalError', ['I''ve considered any ' ...
%         'partition I know and still don''t get the volume of the confinement ' ...
%         'cone. :/ This should not happen.'])
% end